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SUMMARY 


We  compare  a  number  of  advanced  antenna  array  and  time  series  signal 
processing  techniques  which  either  utilise  singular  value  decomposition, 
or  may  be  profitably  interpreted  in  terms  of  it.  Computer  simulations 
of  these  techniques  acting  on  spatio-temporal  data  from  linear,  conformal 
circular  and  focal  plane  arrays  are  presented  and  discussed.  Means  of 
introducing  prior  knowledge  are  described,  and  lead  to  the  suggestion  of 
a  computationally  efficient  processing  approach.  We  have  made  extensive 
use  of  a  sequential  processing  strategy  for  multiple  domain  data  and  show 
by  example  that  this  potentially  allows  both  detection  and  classification 
of  additional  emitters. 
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1.  INTRODUCTION 


The  numerical  technique  known  as  singular  value  decomposition  (SVD)  is 
presented  here  as  a  valuable  tool  (or  many  signal  processing  applications  A 
description  is  given  of  a  number  of  spectral  analysis  techniques  which  either 
utilise  this  fundamental  tool,  or  may  be  usefully  interpreted  in  terms  of  it  By 
comparison  of  these  methods  acting  on  simulated  data  we  demonstrate  the 
usefulness  of  SVD  for  performing  data  conditioning  and.,  as  a  result,  for 
obtaining  improved  signal  estimation. 

Section  2  provides  a  brief  outline  of  SVD  as  a  numerical  tool  for  deriving 
least  mean  squares  solutions  to  ill-conditioned  or  error-sensitive  problems., 
and  shows  how  the  decomposition  gives  a  measure  of  the  contrC.ution  of 
orthonormal  components  to  the  matrix  rank 

In  Section  3.  following  an  explanation  of  the  means  by  which  prior 
information  may  be  introduced  into  the  problem,  a  number  of  methods  ate 
discussed  in  some  detail.  The  methods  described  are  : 

1.  the  Beam  Scan  algorithm,  a  simple  matched  filter  method.. 

2.  the  Pseudo-Inverse  method. 

3.  the  Principal  Components  Bartlett  algorithm.. 

4.  the  Capon  Maximum  Likelihood  method. 

5  the  Thermal  Noise  algorithm. 

6  the  Maximum  Entropy  method. 

7  Multiple  Signal  Classification. 

8  the  Pisarenko  Minimum  Eigenvector  method. 

9.  the  Kumaresan  and  Tufts  Signal  Cancellation  algorithms., 
and  10.  the  Sub-Aperture  technique. 

After  a  discussion  of  the  computational  load  presented  by  the  algorithms., 
we  suggest  possible  means  by  which  a  number  of  them  may  be  accelerated 
This  is  followed  by  the  description  of  a  pre-processing  step  which  makes  use 
of  knowledge  concerning  the  redundancy  in  an  array  of  sensors  in  order  to 
compact  the  data  before  analysis  Consideration  is  also  given  to  the  possibility 
of  utilising  the  output  of  a  number  of  the  algorithms  as  prior  Knowledge  in 
some  form  of  iterative  processing,  or  to  enable  more  detailed  investigation  of 
other  data  domains.  This  leads  us  to  propose  a  computationally  efficient 
sequential  processing  strategy  for  dealing  with  multiole  domain  data.,  which  we 
have  used  extei.  ively  in  our  experiments  However.,  we  show  by  examples 
presented  in  Section  4  that  adoption  of  this  approach  implies  that  the  order  in 
which  domains  of  tne  data  are  investigated  can  have  an  impact  on  the  limits  of 
discrimination. 
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A  wide  variety  of  typical  results  o*  •'omputer  simulations  are  presented  and 
discussed  in  Section  4.  Tnese  demonstrate  the  algorithms  acting  upon 
spatio-temporal  data  generated  by  simplified  models  of  linear,,  conformal 
circular  and  focal  plane  arrays  of  spatial  signal  detectors.  and  are  used  to 
illustrate  many  of  the  points  developed  earlier. 


g.  SINGULAR  VALUE  DECOMPOSITION  (SVD) 


Determination  of  the  singular  value  decomposition  of  a  matrix  is  a  well 
established  numerical  technique  for  evaluating  least-squares  solutions  to 
error-sensitive  or  ill-conditioned  problems  (Lawson  and  Hanson  (1974)..  Golub 
and  Reinsch  (1970)).- 

The  general  matrix..  A.  of  dimensions  (m  x  n) .,  with  m  <  n,,  may  be 
decomposed  as  follows  . 

A  =  U  S  »_H  (2  1) 

where  H  denotes  the  Henmtian  (complex  ccnfugate)  transpose,  and 

yH  v  =  yH  u  =  u  uH  =  .  (2.2) 

G  is  a  Diagonal  matrix  containing  the  non-negative  square-roots  of  the 

si  jenvalues  of  (A  A*"*) .  conventionally  arranged  in  non-i. 'creasing  order  The 

columns  of  the  matrix  V  are  the  m  orthonormalised  eigem-.ctors  corresponding 

to  the  m  largest  eigenvalues  of  (AH  A) .  *  id  the  columns  of  U  are  the 

orthonormalised  eigenvectors  of  (A  AH)  The  columns  of  U  and  of  V  are 

respectively  Known  as  the  left  and  right  hand  singular  '’eotois  of  A.  and  the 

elements  a  .  .  .  ,«r_  of  the  matrix  S  are  called  the  singular  values  of  A  l„  is 

x  m  —  —  “frl 

the  (m  x  m)  identity  matrix 

If  the  matrix..  A  is  square  (m  =  n)  it  is  said  to  be  full  rank  if  the  singular 

values  all  have  value.,  greater  than  zero  M^wever,  if  the  vectors  comprising 

the  rews  or  columns  of  A  are  not  fully  independent.,  then  a  number  of  the 

singular  values  will  tend  towards.,  or  become.,  zero  and  the  matrix  may  no 
longer  be  full  rank  The  rank  of  a  matrix  is  defined  as  being  equal  to  the 
number  of  non-zero  singular  values.,  and  so  for  a  rectangular  matrix  (m  <  n) 
will  always  be  limited  by  the  smaller  dimension.,  m  By  setting  any  very  small 
singuiat  values  to  zero  in  equation  2.1..  we  can  easily  form  a  reduced  rank 
approximation  to  A  Similarly.,  by  inverting  all  non-zero  singular  values,  we 
can  create  a  full  rank  or  reduced  rank  pseudo-inverse  (Penrose  (1955))  of 
the  original  matrix  as  follows: 

A  -  =  V  S'*  yH  (2.  3) 

It  can  be  shown  that  this  procedure  enables  us  to  find  a  least  squares 
solution  to  a  general  linear  problem  (Lawson  and  Hanson  (1974))..  and  may 
therefore  be  considered  as  a  potentially  useful  signal  processing  tool  with 
which  to  average  noisy  data  or  generally  to  solve  numerically  sensitive 
problems  (see.,  for  example.  Varah  (1973).  Cullum  (1980)).  The  theory  may 
also  be  extended  to  deal  with  inverse  problems  involving  both  continuous  and 
discrete  functions.,  and  thus  in  the  solution  of  certain  types  of  integral 
equations  (Hanson  (1971)..  Pike  et  al  (1984)). 
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3.  SPECTRAL  ESTIMATION  ALGORITHMS 


3.1.  SIGNAL  PROCESSING  AND  PftiOR  INFORMATION 


All  signal  processing  relies  on  prior  assumptions  or  information  about  the 
origin  of  the  data.  Presented  simply  with  data  In  the  absence  of  any  other 
contextual  Information,  we  can  say  little.,  if  anything,  about  Its  causation.  In 
order  to  recover  a  transformed  signal  from  received  data.,  we  need  to  know  at 
least  some  of  the  characteristics  of  the  transformation  involved.,  such  as  the 
filter  bandwidth,  and  the  transform  kernel  (spatial  Founer  transform,  tor 
e'  ample) .  We  may  also  make  use  of  information  about  the  general  nature  of 
the  solution  required  -  the  temporal  or  spatial  extent  of  a  signal.,  or  its 
bandwidth  for  example  -  and  the  statistical  form  of  the  noise  in  general., 
signal  processing  should  be  rogarded  as  a  multi-dimensional  problem.,  and  it 
may  be  possible  to  utilise  information  from  one  domain  or  source  of  d3ta  to 
assist  with  the  comprehension  of  other  domains  Such  diverse  information, 
which  is  often  collectively  referred  to  as  prior  knowledge,  not  only  enables  an 
attempt  to  be  made  at  source  reconstruction  if  used  sensibly.,  but  can  also 
give  an  indication  of  the  limiting  factors  Involved  if  it  is  introduced  in  an 
explicit  manner 

The  algorithms  to  be  presented  make  use  of  varying  levels  of  detail  in  the 
possible  range  of  prior  knowledge.,  from  postulation  only  of  the  general  nature 
«_f  the  transform,  to  the  inclusion  of  knowledge  about  the  statistical  behaviour 
of  the  noise.,  and  details  of  the  type  of  solution  required 

3.  g.  MATCHED.  INVERSE  ANO  ORTHOGONAL  FILTERING 

Any  general  linear  system  may  be  represented  as  in  Fig  3  1.  The  input 
vector,  f.  representing  the  source  characteristics.,  is  transformed  by  the 
system  impulse  response..  A.  to  give  a  data  'snapshot'  vector.,  d 

A  f  =  d  (3.  1> 

A  conventional  method  of  signal  processing  attempts  to  transform  tho  data 
into  an  image  of  the  source  by  perform ir  g  the  following  operation  . 

Ah  d  =  f  (3  2) 

This  is  loose  y  known  as  'matched'  or  'correlation'  filtering  :  f  may  be  seen 
as  representing  the  degree  ol  correlation  or  match  between  those  vectors 
which  define  the  system  transform  and  the  data  Generally  defining  the 
columns  of  A  as  representing  the  response  cf  the  system  to  all  possible 
impulses.,  each  row  of  AH  is  a  matched  filter  for  at  least  one  possible  input 
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If  we  take  as  a  specific  example  a  multi-element  linear  at  ray  of  spatial 
signal  detecting  elements,,  then  our  prior  knowledge  may  be  that  A  represents 
a  spatial  Fourier  transform,,  and  AH  Is  the  inverse  Fourier  transformation  This 
type  of  processing  is  only  strictly  valid  if  the  data  vector  is  infinitely  long.,  and 
in  the  case  of  a  Fourier  transform,,  the  effect  of  truncating  the  data  is  to 
create  the  familiar  sin(x)/x  response  to  a  single  frequency  input  signal.  In 
general.;  .‘herefore.  we  might  believe  that  an  exact  inverse  of  A  would  enable 
us  to  recover  f  more  accurately,  thus: 

A-1  d  =  f  (3.3) 

This  is  known  as  inverse  filtering,  and  incorporates  more  precise  prior 
knowledge  of  the  filtering  characteristics  of  the  mapping  A.  This  is  achieved 
through  the  implied  constraint  on  the  extent  of  f  which  is  required  for  the 
formation  of  the  inverse  For  a  sample  of  data  from  an  infinitely  long  line 
array.,  this  becomes  equivalent  once  again  to  the  Inverse  Fourier  transform,, 
since  A  AH  =  Ijj,  in  this  case. 

Orthogonal  filtering  involves  comparing  the  system  response  matrix  with  a 
vector  which  is  orthogonal  to  the  data.  In  the  case  of  the  spatial  sensor  array., 
the  result  is  a  vector  with  nulls  In  the  directions  of  decorrelated  sources.  This 
vector  is  conventionally  plotted  inverted  on  a  logarithmic  scale.,  to  give  peaks 
at  the  appropriate  angles.  An  example  of  such  a  filter  would  be  an  adaptive 
beamformer  network  (Applebaum  (1976))  which  is  designed  to  filter  out  the 
spatial  frequencies  in  the  data  which  correspond  to  sources  of  interference 

A  number  of  the  algorithms  to  be  described  make  use  of  all  the  three 
principles  outlined  above,  utilising  the  results  from  the  orthogonal  filter  to 
provide  support  for  the  design  of  an  Inverse  filter  which  is  accurately 
constructed  from  vectors  matched  to  the  signals  detected.  However.,  before 
describing  specific  examoles  of  algorithms,,  the  system  response  matrix,  A.  will 
be  discussed  in  greater  detail.  For  a  spatial  array  of  sensors,  this  is 
commonlv  known  as  the  array  manifold  as  a  result  of  its  multi-dimensional 
vector  space  interpretation  (Schmidt  (1981)). 

3.3.  THE  ARRAY  MANIFOLD 

3,3.1.  CALIBRATION 

The  array  manifold  is  a  multi-dimensional  characterisation  of  the  impulse 
response  of  the  particular  configuration  of  sensors  under  consideration.  We 
show  first  how  this  characterisation  may  bo  carried  out  for  the  example  of  a 
linear  phased  array  antenna.  The  method  Is  easily  applied  to  alternative 
conformations. 
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As  depicted  by  Fig  3. 2.  a  signal  source  in  the  far  field  will  result  In  a 
plane  w«rve  being  received  at  a  particular  angle,  (rr/2  ±  e) .  to  the  fdC9  of  the 
array.  Far  example,  a  source  along  a  direction  normal  to  the  plane  of  the 
array  (e  =  0)  will  result  In  a  wavefront  parallel  to  the  array.  The  simplifying 
assumption  of  uniform  amplitude  and  equal  phase  weighting  applied  to  the 
sensing  elements  of  the  array  will  then  produce  a  data  vector  sampled  at  an 
Instant  in  time  whose  terms  are  all  equal.  Similarly.,  a  source  mcated  at  a 
positive  (negative)  angle  to  the  normal  will  result  in  positive  'negative) 
frequency  spatial  pattern  i.  the  data  vector,  in  the  simplified  case  this  pattern 
is  a  discrete  complex  exponential  sampled  at  the  element  spacings  and 
truncated  by  the  spatial  limits  of  the  array. 

Assuming  that  sources  are  in  the  far  field  and  approximate  to  point 
emitters,  complex  data  vectors  may  be  recorded  for  all  such  possible  source 
locations  within  the  field  of  view  of  the  antenna.  These  vectors  can  then  be 
collected  together  as  the  columns.  a(e).  of  a  matrix.  A.  which  forms  a 
'reference  library'  of  signals  expected  for  each  incident  direction.  This  library 
is  a  discrete  space  equivalent  of  the  array  manifold  described  by  Schmidt 
(1981).  Since  each  of  the  vectors  which  define  the  manifold  may  also  be  used 
to  control  the  phases  of  each  antenna  element  in  a  beamforming  mode.,  they 
are  often  referred  to  as  'steering  vectors'. 

3-3.2.  AM8IQUITIES 

Clearly,  ambiguities  may  ex-st  in  the  interpretation  of  the  steering  vectors 
(Schmidt  (1979)).  Even  assuming  that  spatial  sampling  rate  criteria  have  been 
met.  in  the  instance  of  a  linear  array  of  receiving  elements  as  described.,  the 
array  manifold  will  still  contain  one  ambiguity  .  every  column  In  the  library  may 
be  interpreted  either  as  the  response  to  a  source  ahead  of  the  array.,  or  to 
one  behind  the  array  This  ambiguity  may  be  resolvable  by  the  addition  of 
furiher  Information  about  the  array.,  such  as  that  signals  are  only  to  be 
expected  from  sources  ahead  of  the  array.,  or  that  the  fie«o  of  view  is  stil1 
further  restricted  by  the  directional  properties  of  each  element  for  example. 
The  field  of  view  may  be  specified  In  a  general  way  by  staling  its  extent  in 
terms  of  a  number  of  Rayleigh  beamwidths  (Cox  (1973))  defined  primarily  Dy 
the  overall  3perture  of  the  antenna.  This  is  analogous  to  the  information 
transmission  capacity  of  the  system  (Shannon  (1948)),.  and  is  described  by  its 
Shannon  number. 

Having  demonstrated  the  way  In  which  we  may  construct  the  array 
manifold,  an  example  of  the  real  (cosine)  component  of  such  a  manifold  is 
given  in  Fig  3  3.  The  array  is  a  co-linear  phased  array  of  twenty  elements 
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Fig  3  3  The  real  (cosine)  component  of  a  typical  array  manifold  Each 
column  represents  tha  array  response  to  a  source  at  ona  of  60  angles  in 
spaca  This  particular  20  element  array  has  a  Shannon  number  of  3  so  that 
tha  limits  to  the  angle  of  view  are  represented  by  waveforms  of  1  5  cycles 
lii«t  ambiguity  between  left  and  right  is  resolved  by  the  imaginary  component 
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with  uniform  antenna  weights  and  spacing.,  and  having  a  Shannon  number  of 
six.  The  field  of  view,  or  object  domain.,  is  quantised  into  sixty  positions 
across  the  six  beamwidths.  A  measure  of  the  angle  of  arrival  of  signals  is  the 
number  of  spatial  cycles  created  across  the  face  of  the  array  by  the  source 
For  this  case  the  angular  limits  of  the  field  of  view  thus  correspond  to  *  3 
cycles. 

3.3.3.  THE  MANIFOLD  AS  PRIOR  KNOWLEDGE 

Conventionally..  Interpretation  of  received  data  from  a  line  array 
corresponds  directly  to  performing  a  discrete  Fourier  transform.  For  a  situation 
in  which  the  array  response  does  not  approximate  closely  to  the  idealised 
linear  example  given  above.  Fourier  techniques  are  not  directly  applicable  The 
methods  of  analysis  to  be  described  are  designed  to  take  account  of  any  such 
perturbations  if  they  can  be  modelled  or  measured. 

It  can  clearly  be  seen  that  the  array  manifold  contains  much  prior 
knowledge  useful  for  interpreting  the  received  data.  Firstly,  it  holds  the  system 
response  characteristics  for  a  large  number  of  directions  of  arrival.,  and 
secondly  it  may  also  implicitly  contain  knowledge  about  the  parameters 
associated  with  particular  signals.,  or  the  angular  bounds  between  which  signals 
are  to  be  expected. 

3. 4.  THE  SAMPLE  COVARIANCE  MATRIX 

3.4.1.  UNCORRELATED  SOURCES 

If  we  have  a  sequence  of  data  snapshots  of  a  source  scenario  consisting 
of  uncorrelated  point  emitters  It  should  be  possible  to  make  use  of  the 
additional  Information  that  the  temporal  modulation  in  the  data  provides  It  can 
be  shown  that  this  extra  information  may  be  optimally  recovered  by  examining 
the  eigenvalue  spectrum  of  the  data  covariance  matrix  (van  Trees  (1968)) 

Assuming  that  signals  and  noise  are  uncorrelated.,  and  that  the  noise  is 
zero  mean  Gaussian  at  each  antenna  element,  we  represent  the  process  of 
detection  as  : 

O  =  A  F  ♦  N  (3  4) 

where  D  is  the  matrix  of  snapshots.,  the  columns  of  £  contain  the  unknown 
target  parameters  as  they  vary  with  time.,  and  N  is  the  record  of  noise 
samples. 

For  uncorreiated  signals  with  Gaussian  properties.,  knowledge  of  the 
correlation  function  of  the  process  Is  sufficient  to  define  the  required 
spectrum.  However,,  in  general  we  do  not  have  sumcient  samples  to  form  an 
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unbiased  estimate  of  this  function,  and  so  we  compute  the  data  covariance 
matrix,  which.,  in  the  large  sample  limit  tends  towards  the  form  of  a  ToepliU 
correlation  matrix  (Makhoul  (1975)). 

The  (m  x  m)  spatial  covariance  matrix  of  the  data  may  be  estimated  in 
practice  by  a  variety  of  techniques,  such  as  simple  block  averaging.,  sliding 
window  averaging  (see  section  3.7)..  time  decay  averaging,  forward-backward 
averaging.,  etc  (Gabriel  (1984))..  Since  these  methods  may  all  be  Interpreted 
as  differences  in  the  construction  of  the  data  matrix.,  they  can  be  generally 
defined  by  the  sample  covariance  estimate.,  R.  given  by  the  equation,. 

R=DDH  =  AHAH  +  g  =  ttf  +  g  (3  5) 

where  H  is  F  FH.  the  'signal-in-space'  covariance  matrix..  W  is  the  signal 
covariance  matrix,  and  Q  is  the  noise  covariance  matrix.  Over  a  long  period  of 
observation  the  uncorrelated  nature  of  the  noise  will  result  in  g  becoming  a 
diagonal  matrix,  which  will  simply  add  to  the  eigenvalues  of  the  matrix  W.  If 
the  noise  is  white.,  equation  3.  5  may  be  written. 

R  =  W  +  cr  *1  .  (3.6) 

—  O  “ 

where  I  is  the  Identity  matrix,  and  is  tho  noise  power  which  will  equal  the 
minimum  eigenvalue.  Amjn<  of  the  generalised  eigen  equation. 

B  ♦,  =  *,  9  ♦,  1  <  K  m  ..  (3. 7) 

3.  4.  g-  THE  SIGNAL  AND  NOISE  SUBSPACES 

The  signal  covariance  matrix  may  be  considered  as  a  sum  of  outer 
(dyadic)  products  of  the  columns  of  the  A  matrix,  weighted  according  to  the 
power  of  the  corresponding  signal  associated  with  each  direction.  Thus.,  since 
the  rank  of  each  outer  product  is  one.  the  rank  of  the  total  signal  correlation 
matrix  is.  in  the  absence  of  noise,  and  given  a  sufficient  number  of 
independent  data  snapshots,  equal  to  the  number  of  spatially  and  temporally 
uncorrelated  sources.  Ns.  for  Ng  <  m.  In  the  presence  of  noise,  an 
evaluation  of  the  pseudo  rank  of  R  must  be  made  in  order  to  form  an  estimate 
of  the  minimum  number  of  sources  which  are  needed  for  the  model. 

In  general,  if  we  have  a  good  measurement  of  the  second  order  noise 
statistics  in  the  form  of  the  matrix  Q.  we  may  solve  for  the  eigenvalues  of  R  in 
the  metric  of  the  noise  (equation  3.7).-  Cox  (1973)  shows  that  this  is 
equivalent  to  ensuring  that  emphasis  is  given  to  those  components  containing 
least  noise  so  that  the  choice  of  pseudo  rank  for  P  becomes  less  subjective. 
The  eigenvectors  of  R  associated  with  the  largest  eigenvalues  are  known  as  the 
'signal  subspace*  basis  vectors,  and  those  associated  with  the  smaller 
eigenvalues  are  referred  to  as  the  'noise  subspace'  basis  vectors.  Because  of 
the  mutual  orthogonality  of  singular  vectors  (the  eigenvectors  of  the  covariance 
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matrix  are  the  singular  vectors  of  the  data)  the  two  subspaces  are  also 
mutually  orthogonal.  We  therefore  appear  to  have  a  means  by  which  the 
signals  and  noise  may  be  largely  separated  if  required.  It  is  important  to 
remember,  however,,  that  the  signal  and  noiso  as  such  are  not  necessarily 
completely  orthogonalised..  as  the  subspace  terminology  might  be  taken  to 
imply.  Given  a  limited  sample  of  data.,  it  may  be  the  case  that  the  signal 
subspace  eigenvectors  still  contain  a  significant  component  of  noise,,  and  vice 
versa 

3.4.3,  CORRELATED  SOURCES 

Implicit  in  the  above  discussion  is  the  assumption  that  the  autocorrelation 
function  is  stationary,,  which  is  the  same  as  requiring  that  the  covariance 
matrix  tends  to  the  Toeplitz  form  in  the  large  sample  limit.  This  is  only  the 
case  if  it  is  true  both  that  the  sources  ate  stationary  in  space  for  the  duration 
of  the  observation  and  that  they  appear  to  be  decorrelated  within  that  interval. 

Correlated  sources  produce  fields  in  space  which  are  non-stationary  (White 
(1979))  and  the  data  matrix  which  is  formed  does  not  have  a  pioportional 
relationship  between  rank  and  number  of  emitters.  Thus  in  general,  processing 
intended  for  the  estimation  of  stationary  signals  does  not  produce  useful 
results,  unless  the  data  may  be  pre-processed  In  such  a  way  as  to  derive  an 
approximation  to  a  stationary  correlation  function.  A  possible  means  of 
achieving  this  is  described  in  section  3.7. 

3.4.4.  THE  MINIMUM  DESCRIPTION  LENGTH 

Numerous  methods  have  been  devised  for  automatically  evaluating  the 
dimension  of  the  signal  subspace  by  examination  of  the  eigenvalues  of  R 
These  variously  make  use  of  curve  fitting  techniques  (Gabriel  (1984)).. 
maximum  likelihood  estimates  of  the  model  order  (Akaike  (1974)).  and 
Information  theoretic  criteria  (Rissanen  (1978)).  We  have  chosen  to  implement 
the  minimum  description  length  (MOL)  criterion  based  on  the  work  of  Rissanen 
(1978)  and  described  by  Wax  and  Kailath  (1983).  This  criterion  forms  the 
estimate  by  calculating  the  function 

MDL(k)  =  -  Nd  log{  n  a^/ll/fm-k)  £  crj*)m-k)  ♦  ’*k(2m-k) log( N0>  (3  8) 

where  Np  is  the  number  c!  observations  or  snapshots.,  for  k  =  0  to  m-1..  and 
finding  the  value  of  k  for  which  it  is  minimised.  The  first  term  of  this  equation 
is  th6  log  of  the  maximum  likelihood  estimates  of  the  model  parameters  (van 
Trees  (1968))  and  the  second  Is  a  bias  correction  term.  This  criterion  has 
been  used  in  our  simulations  and  has  consistently  given  estimates  of  subspace 
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Fig  3  4a  A  typical  amplitude  function  consisting  of  a  turn  of  two  tine 
functions  demonstrating  'classical*  resolution. 


Fig  3.4b 
amplitude 


A  typical  function  consisting  of  two  sine  functions  of  unequal 
in  a  classical  sense,  resolution  has  not  baan  achieved. 


dimension  which  agree  with  those  which  would  have  been  selected  by  a  manual 
operator, 

3.  5.  RESOLUTION  AND  DISCRIMINATION 

Cox  (1973)  defines  'resolution'  in  the  classical  sense  simpiy  as  recognition 
that  an  observed  effect  Is  due  to  two  separate  sources  instead  of  one. 

practice.,  perhaps  the  best  we  can  say  is  that  the  observation  is  the  result  of 

a  minimum  of  two  sources,,  since  the  effect  of  additional  sources  may  be  lost 
In  the  noise.  However,  the  definition  elaborated  by  Cox  Is  based  on  the 
assumption  that  resolution  may  be  measured  by  the  ability  of  a  processor  to 
produce  distinct  peaks  In  its  response  corresponding  to  each  independent 
source  of  signal  (Fig  3  4a).  This  Immediately  implies  a  direct  'operator' 
interpretation  of  the  spectrum  in  terms  of  point  targets  corresponding  to  such 
peaks. 

In  the  example  of  an  ideal  regular  linear  array,  we  know  that  this  is  an 
incorrect  interpretation  of  the  end  product  of  Fourier  processing.  Prior 
knowledge  of  the  transform  kernel  In  this  case  tells  us  that  the  correct 
interpretation  is  as  a  sum  of  sine  functions  of  unknown  amplitudes  and 

positions,  in  this  sense  we  are  able  to  say  that  the  Fourier  beamformer  has 

resolved  when  It  is  apparent  that  the  response  differs  from  the  single  sine 

function  in  some  way  such  as  illustrated  in  Fig  3. 4b.  Clearly  this  does  not 

lead  to  a  very  useful  definition  of  resolution  in  practice  since  it  gives  us  no 
direct  measure  of  relative  target  position.,  other  than  that  the  sources  are 
placed  within  the  conventional  Rayleigh  limit.  We  therefore  prefer  to  use  the 
term  'discrimination'  to  describe  the  ability  of  an  algorithm  to  produce  distinct 
peaks  or  nulls  for  each  source  in  its  output. 

3.6.  ALGORITHMS 

3.6.1.  THE  BEAM  SCAN  ALGORITHM  (BSA) 

The  beam  scan  algorithm  is  otherwise  loosely  known  as  the  correlation  or 
matched  filter.,  since  it  is  composed  of  a  number  of  parallel  filters  (the  rows 

of  Ah>  ..  each  matched  to  a  specific  possible  Input.  We  refer  to  It  as  a  beam 

scan  algorithm  because  the  signal  processing  is  analogous  to  scanning  the 
beam  of  a  radar  array  through  the  data.  The  data  is  compared  with  each  of 
the  beam  steering  vectors  contained  in  the  array  manifold  in  turn.,  and  the 
resulting  product  forms  the  spatial  power  estimate.  If  we  have  more  than  one 
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data  vector  (multiple  snapshots).,  the  sample  covariance  estimate..  R.  may  be 
calculated.. 

R  =  D  DH  .  (3. 9) 

and  the  power  estimate  may  be  written  as. 

P__..(e)  =  aH(e>  Rate)  (3.10) 

where  e  is  the  angle  of  arrival  of  the  signal,  and  a(e)  Is  the  corresponding 
steering  vector. 

The  ability  of  the  BSA  to  discriminate  multiple  spatially  distributed  point 
emitters  by  producing  Independent  peaks  In  Pqsa(8>  Is  limited  for  short  data 

records  because  of  the  resulting  broad  'beam'  shape  and  high  sidelobes..  but 

o 

nevertheless  provides  the  best  possible  (maximum  likelihood)  linear  estimator 
of  the  direction  of  a  single  emitting  source  from  data  containing  additive  white 
Gaussian  noise  (van  Trees  (1968)). 


3.6.2.  THE  PSEUDO-INVERSE  METHOD  (PIM) 

As  stated  in  section  3.2..  an  improved  processor  for  a  given  system  should 
be  an  exact  inverse  of  the  impulse  response  matrix.  A  in  order  to  recover  f 
thus:. 

A‘*d  =  f  (3.  11) 

In  the  case  of  »he  linear  detector  array.,  because  the  vectors  contained  in 
the  library  matrix  described  in  section  3. 3  consist  of  truncated  complex 
exponentials,  they  are  not  necessarily  fully  Independent.  This  results  in  the 
inversion  of  A  being  an  ill-conditioned  problem  since  a  number  of  the  singular 
values  become  very  small.  The  operation  described  by  equation  3.11  implies 
the  multiplication  of  components  within  the  data  by  the  inverted  singular 
values.  For  the  transform  kernels  of  interest  to  us  here  this  would  possibly 
lead  to  the  amplification  of  any  noise  contained  in  the  data  and  of 
uncertainties  in  the  characterisation  of  the  array  (Cullum  (1980)).  The  array 
manifold  must  therefore  be  approximated  In  some  suitable  way  The 
computational  problem  Is  further  exacerbated  by  the  generalisation  to  a 
rectangular  A  matrix,  which  represents  a  limited  number  of  sample  points  in 
relation  to  an  infinite  number  of  possible  source  locations.,  and  may  be  thought 
of  as  resulting  in  an  underdetermined  system  of  equations 

One  means  of  overcoming  both  of  these  problems  is  to  recombine  only  a 
small  number  of  singular  values  and  vectors  from  the  decomposition  of  A  in 
order  to  provide  a  least-squares  solution  to  *he  problem  (Lawson  and  Hanson 
(1974)).  A  suitable  upper  limit  on  the  number  of  singular  vectors  used  could 
occur  when  the  calculated  singulai  values  fall  to  the  level  of  errors  introduced 
by  their  computation.  In  the  example  of  the  twenty  element  array  used  earlier.. 
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with  Its  Shannon  number  of  six.  truncating  the  singular  value  series  when  ci/o1_ 
fails  below  say  10-s  leaves  us  with  only  eleven  'significant'  singular  vectors 
from  the  initial  set  of  twenty.  The  array  may  be  said  to  have  only  eleven 
significant  degrees  of  freedom  (Toraldo  di  Francia  (1969)).  Truncation  ot  the 
singular  value  series  now  allows  us  to  find  a  pseudo-inverse  in  the  general 
case  of  a  rectangular  A  matrix 

Using  the  singular  value  decomposition  of  A. 

A  =  U  S  V1"1  (3  12) 

we  may  re-write  equation 

A  f  -  d  (3  13) 

in  the  form 

S  VH  f  =  UH  d  =  y  .  (3  14> 

by  decomposing  the  data  onto  the  orthonormai  basis  provided  by  the  columns 
of  U  Finally  the  result  is  obtained  from 

f  =  V  £  <y|/c,)e,  =  V  S-1y  (3  15) 

i*  i 

where  the  y,  are  components  of  y.  and  (e()  are  the  columns  of  an  identity 
matrix  Thus.. 

f  =  v  S-1  UH  d  (3. 16) 

Clearly,  if  we  are  finding  a  rank-reduced  pseudo-inverse  of  A  by  setting  some 
of  the  cj  to  zero.,  we  must  assume  that  the  y,  decrease  faster  than  do  the  a, 
with  increasing  I.  If  this  Is  the  case  we  may  then  truncate  the  a,  when  they 
fall  to  the  level  of  the  noise  in  the  system,  since  this  implies  that  eigenvector 
components  of  the  data  have  been  attenuated  by  the  transformation  to  such  a 
degree  that  they  may  no  longer  be  recovered. 

For  multiple  snapshots  of  data.,  the  power  estimate  as  a  function  of  angle 
of  arrival  may  be  written  in  the  same  form  as  for  the  BSA., 

Pp(M(e)  =  a-1(e>  R  a'(e)  (3  17) 

where  a-1(e>  is  a  row  of  the  matrix  A"1.,  and  a’(e>  is  a  column  of  the 
rank-reduced  system  transform 

The  'beamwidth'  of  this  processor  varies  according  to  the  number  o* 
singular  values  included  in  the  pseudo-inverse  calculation,  and  under 
conditions  of  low  signal  to  noise  ratio.,  with  appropriate  truncation  of  the  o,. 
broadens  until  it  approaches  that  of  the  BSA,  Very  high  S/N  ratios  have  to  be 
sustained  in  order  to  achieve  beamwidth  reductions  of  up  to  a  factor  of  3  or 
4.  and  the  algorithm  is  much  more  sensitive  to  the  validity  of  the  assumed 
spatial  extent  of  the  object  domain  than  the  BSA.,  as  will  be  shown  by 
examples  presented  in  Section  4 
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3.6.3.  THE  PRINCIPAL  COMPONENTS  BARTLETT  METHOD  ( PCS) 


This  signal  matching  method  estimates  the  spatial  power  spectrum  by  a 
similar  approach  to  the  SSA  It  differs  from  BSA  in  using  a  rank-reduced  oata 
matrix  as  the  object  to  be  scanned  by  the  beam 

In  section  3  4  it  was  shown  that  the  sample  covariance  matrix  cou'd  be 
diagonallsed  in  older  to  define  a  signa>  subspace,  with  basis  vectors 
corresponding  to  those  eigenvectors  with  large  eigenvalues,  and  an  orthogonal 
noise  subspace  defined  by  those  eigenvectors  associated  with  sma'l 
eigenvalues  The  PCB  algorithm  involves  the  formation  of  an  estimate  of  R 
which  incudes  only  the  signal  space,  or  principal  eigenvectors,  and  sets  the 
corresponding  eigenvalues  to  unity, 

Rs  =  ’jf'*,  «-(H  (3  18) 

Inherent  in  this  description  is  a  strong  interpretation  of  the  eigenvalue 
spectrum  VV6  are  asserting  that  the  large  eigenvalues  do  correspond  to 
signal,  and  that  the  rest  do  correspond  to  noise  as  far  as  subsequent 
processing  is  concerned  Particularly  in  the  case  of  under-estimation  of  the 
Signal  subspace  dimension  incorrect  partitioning  of  the  spaces  can  have 
significant  influence  on  the  form  of  the  spectrum  obtained 

m  the  formation  of  Rg  the  PCB  algorithm  has  thrown  away  the  estimate  it 
had  about  the  relative  signal  powers  in  the  form  of  the  discarded  eigenvalues 
This  may  even  be  considered  a  valid  step,  particularly  in  the  case  of  short 
data  records,  when  tne  estimate  may  have  been  a  poor  one,  and  is  in  fact 
the  means  by  which  »ne  algorithm  ensures  discrimination  of  low  amplitude 
emitters  m  proximity  to  sources  of  higher  power  By  fitting  functions  of  equal 
magnitude  to  the  data,  the  algorithm  circumvents  the  difficulty  in  the  definition 
of  resolution  described  by  Cox  (1973)  which  arises  under  such  circumstances 
We  may  now  simply  form  an  evaluation  of  the  directions  of  arrival  of 
signals, 

PpCB(0)  =  aH(e)  Rg  ate)  (3.19) 

If  we  now  wish  to  interpret  the  resulting  spectrum  in  terms  of  point  targets,  it 
is  necessary  to  locate  the  Np(<  peaks  of  the  generally  very  broad  estimate. 
Having  achieved  this,  the  steering  vectors  corresponding  to  such  peaks  are 
taken  as  indicating  the  directions  of  the  point  emitters,  and  may  be  stored  in 
a  full  rank  matrix.  B  Forming  a  pseudo-inverse  of  B  and  multiplying  by  the 
data  matrix, 

T  =  <Bh  B)'a3h  D  =  B-1  O  (3  20) 

produces  a  matrix  containing  Np*  rows,  each  of  which  is  a  time  series  of  data 
associated  with  a  particular  direction  in  space,  information  from  the  other 
directions  being  cancelled  by  the  r.lternative  steering  vectors  in  B  As 
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described  In  section  3.7.,  further  processing  of  T  can  provide  more  detailed 
information  about  the  target  scenario,  such  as  the  number  of  emitters  at  each 
direction.  Target  cross  and  auto-powers  can  be  derived  by  examining  the 
elements  of  the  reduced  target  data  covariance  matrix.,  CT  TH)  The  powers 
determined  in  this  fashion  are  the  total  powers  corresponding  to  each  direction 
of  arrival  and  may  be  of  use  in  discriminating  against  spurious  peaks  in  the 
response  which  are  the  result  of  noise. 

3.6.4.  THE  CAPON  MAXIMUM  LIKELIHOOD  METHOD  (MLM) 

The  maximum  likelihood  method  (Capon  et  al  (1967))  may  be  interpreted 
as  an  orthogonal  filtering  technique.,  since  it  utilises  the  inverse  of  the  sample 
covariance  matrix  to  produce  s  vector  with  nulls  in  the  di-ection  of  sources, 
and  with  null  depths  corresponding  to  estimates  of  signal  strengths.  The 
spatial  spectral  estimate  is  then  derived  from  the  inverse  of  the  component 
magnitudes  of  the  vector  containing  the  nulls. 

The  approach  to  this  algorithm  is  common  to  all  the  orthogonal  methods  to 
be  discussed.  An  attempt  is  made  to  minimise  the  mean  square  array  output 
power  under  certain  constraints  It  is  the  constraints  which  differentiate  the 
various  algorithms  In  the  case  of  MLM  the  constraint  is  that  the  array  should 
have  unit  gam  in  a  particular  direction.,  whilst  simultaneously  performing  the 
power  minimisation.  This  has  the  effect  of  optimally  rejectirg  power  trom 
directions  other  than  that  which  defines  the  constraint  (Cox  (  1973)). 

For  the  case  of  Gaussian  noise,  the  average  power  output  from  the  array 
for  a  plane  wave  arriving  from  an  angle  e  is  given  by 

PMLM(e)  =  w(e)H  B  w(e)  •  (3  21) 

where  w(e)  is  the  vector  of  array  weights.  The  unity  gain  constraint  requires 
that  w(e)  should  satisfy 

w(e)H  a(e>  =  l  (3.22) 

The  optimum  w  vector  is  the  Wiener  weight  vf  tor.,  which  maximises  the  output 
signal  to  interference  ratio. 

w(e)  =  k  FT1  a(e)  .,  (3  23) 

where  k  Is  a  complex  number.  Solving  for  k  in  equation  3  23..  we  find  that  k 
and  Pmlm  are  (be  same.  Substituting  for  k  and  solving  for  w(e)  produces  the 
MLM  power  estimate. 

P„,  „(e)  =  (  a(e)H  FT1  a(»>  )_1  (3  24) 

Interpreting  this  result  In  terms  of  an  eigenvector  decomposition  ot  R.  the 
matrix  inverse  is  calculated  using  ail  ot  the  eigenvector  and  eigenvalue 
information.  There  is  no  attempt  to  Impose  a  signal  subspace  interpretation 
onto  the  decomposition  The  result  may  be  understood  er.  an  attempt  to  create 
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a  filter  which  is  maximally  orthogonal  to  the  signal.,  since  noise  eigenvectors 
will  receive  greater  weight  in  the  calculation  of  R"1.  without  •.  ctually  having  to 
decide  where  to  partition  the  signal  and  noise  spaces. 

It  should  be  noted  that  the  estimate  in  equation  3.  24  is  the  true  maximum 
likelihood  estimate  only  for  the  assumption  of  a  single  plane  wave.  The 
solution  fcr  multiple  plane  waves  becomes  progressively  more  complicated.,  as 
discussed  by  Evans  et  al  (1982).  Statistical  estimation  theory  shows  that  the 
true  maximum  likelihood  estimate  is.  the  optimal  estimate  for  the  Gaussian 
signal  detection  problem.,  and  enables  discrimination  which  approaches  the 
theoretical  limd  for  sucn  unbiased  estimators,  -refined  by  <lmits  «••...  as  the 
Cramer-'tao  bound  on  the  variance  of  the  estimate  (van  Trees  (1°'8)) 

The  MLM  estimate  has  very  low  'sidelobe'  characteristics  foi  the  case  of 
incoherent  signals,  since  tr(a)  places  nulls  at  all  directions,  unless  e 
corresponds  to  the  direction  of  a-rivai  of  a  plane  wave,  when  the  unity  gain 
constraint  prevents  a  null  be'ng  produced.  The  magnitude  of  the  response  then 
gives  an  estimate  of  tho  power  of  the  amitter. 

3.6.5.  THE  THERMAL  NOISE  ALGORITHM  <TNA) 

in  attempt  to  devise  an  algorithm  which  had  the  desirable  characteristics  of 
low  'sidelobe'  response  coupled  with  sharp  'resolution'  peaks.  Gabriel  (1380) 
derived  his  thermal  noise  algorithm. 

Taking  the  dot  product  of  the  optimum  weight  vector  derived  above  with 
itself,  and  defining  this  as  the  adapted  thermal  noise  power  N0. 

N&  =  A  w(e)H.w(e)  (3.25) 

where  x  Is  a  qc'ascent  noise  power  level  constant,  we  der‘ve  the  new  spectral 
estimator 

PTIMA  =  (  a(e)H  B~z  «(0)  >_1  (3.26) 

from  the  reciprocal  of  N0.  The  peaks  produced  by  this  algorithm  are 
proportional  to  the  square  of  the  source  power  levels. 

if  we  examine  this  method  as  we  did  for  MLM,  by  looking  at  its  use  of  the 
eigenfunctions  of  R.  w9  see  that  the  spatial  estimate  may  be  written 

PTNA<e)=  (  ®(e,H  *  A-*  a(e)  )  (3.27) 

where  ♦  Is  the  matrix  of  eigenvectors  of  R.  and  A  Is  the  diagonal  matrix  of  its 
eigenvalues.  The  TNA  may  therefore  be  interpreted  as  placing  a  stronger 
weighting  on  the  eigenfunctions  which  are  most  likely  to  belong  to  a  noise 
subspace.  As  with  MLM.  the  decision  as  to  the  osrtltionlng  of  the  subspaces 
has  been  avoided,  but  in  this  case  the  implication  that  the  small  eigenvalues 
correspond  to  tne  noise  is  reinforced  by  the  squaring  of  the  inverted 
eigenvalues. 
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Clearly  the  Idea  of  manipu'ating  the  effective  eigenvalue  spectrum  in  this 
way  may  be  extended  to  enable  the  processor  to  impose  even  stronger 
assumptions  about  the  nature  of  the  data  on  the  result.  The  MUSIC  and  PCS 
algorithms  are  examples  of  a  step  function  weighting  being  applied  to  the 
eigenvalues. 

3.6.6.  THE  MAXIMUM  ENTROPY  METHOD  (MEM) 

The  maximum  entropy  method  was  developed  as  a  means  of  extrapolating  a 
limited  estimate  of  an  autocorrelation  function  in  such  a  way  that  the 
appropriate  probability  density  function  has  Its  entropy  maximised  at  each 
stage.  Its  intention  is  to  form  a  'maximally  noncommittal'  estimate  (Jaynes 
(1957)).  In  common  with  MLM,  TNA.  MUSIC.  KTSC.  and  many  other  similar 
algorithms.,  however.,  MEM  is  equivalent  to  a  constrained  least  squares  fitting 
of  an  all-pole  model  to  the  available  data  (van  den  Bos  (1971))  and  may  be 
outlined  without  recourse  to  information  theoretic  concepts  for  the  simple  case 
of  the  regular  linear  array. 

For  brief  data  records,  the  usual  unbiased  estimator  is  tne  data 
covariance  matrix,,  and  we  obtain  the  required  filter  coefficients  by  solving  the 
equation 

R  W  =  [  e  0  0  .  .  .  0  ]T  =  e  (3  28) 

where  e  is  the  SO-called  prediction  error  (Makhoul  (1975)),  with  the  weight 
vector, 

w  =  (  wA  Wj  w,  .  .  wm  )T  (3  29) 

We  obtain  the  solution  for  w  of  the  form 

w  =  R'1  e  (3  30 

and  scaling  w  and  e  appropriately  we  may  choose  e  to  be  the  first  column  of 
the  Identity  matrix.  Thus  the  filter  w  is  equal  to  the  first  column  of  the  inverted 
covariance  matrix  (Nuttall  (1976)).- 

interpreting  this  filter  in  terms  of  eigenfunctions,  as  we  did  with  MLM,  it 
can  be  seen  that  we  again  have  derived  a  function  which  gives  greatest  weight 
to  the  noise  components  of  the  data  without  the  requirement  for  any  additional 
information  or  assumptions  about  the  subspace  dimensions  If  we  compare  this 
filter  against  the  library  of  'expected'  waveforms  in  the  A  matrix,  nulls  should 
be  produced  In  the  result  at  directions  corresponding  to  the  signal  directions. 
Thus  our  spectral  estimate  may  be  written, 

Pmem<®>  =  <  a(e>H  w  wH  a(e)  )~x  (3.31) 

it  Is  worth  noting  the  similarity  In  form  between  this  estimator  and  Ptna-  in 
which  R-^e  Is  replaced  by  R-1.  In  fact  the  TNA  result  may  be  interpreted  as 
an  average  of  MEM  estimates  derived  from  alternative  columns  of  the  inverse 
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covariance  estimate.  Burg  (1972)  has  also  demonstrated  a  link  between  MLM 
and  an  average  of  MEM  estimates  resulting  from  different  filter  lengths,  under 
certain  circumstances.  Note  also  that  the  form  of  equation  3.31  is  only  strictly 
the  maximum  entropy  solution  for  a  regularly  sampled  process  and  that  the 
magnitude  of  Pmem(0)  does  not  necessarily  relate  to  the  received  powers 

3.6.7.  THE  MULTIPLE  SIGNAL  CLASSIFICATION  ALGORITHM  (MUSIC) 

Schmidt  (1979)  has  described  a  signal  subspace  algorithm  which  he  calls 
MUSIC  (Multiple  Signal  Classification).  This  is  fully  described  by  him  for  the 
generalised  case  where  the  data  covariance  matrix  is  analysed  in  the  metric  of 
the  previously  estimated  noise  (equation  3.7).  Clearly  the  noise  covariarcs 
matrix  must  be  measured  In  the  absence  of  other  signals  and  then  assumed  to 
be  stationary  for  a  further  period  whilst  the  generalised  eigenvalue  problem  is 
being  solved  for  each  block  of  data.  Alternatively..  In  order  to  reduce  the 
required  computation,  the  noise  may  be  assumed  Gaussian.,  and  an  estimate 
of  rank  made  from  the  decomposition  of  R  in  an  Identity  metric. 

Having  obtained  the  singular  vectors  of  the  data,  and  formed  an  estimate 
of  the  number  of  identifiable  independent  sources  present  from  the  matrix  rank 
estimate,  the  mutually  orthogonal  signal  and  noise  subspaces  can  be  defined. 
Defining  to  be  the  matrix  of  noise  subspace  basis  vectors,  the  directions  of 

arrival  are  estimated  by  computing  the  Euclidean  distance,  d,  between  the 
vectors  of  the  array  manifold  and  the  noise  subspace  and  then  plotting  1/d2 
against  angle  of  arrival,  e  : 

PMg(e)  =  (a(e)H  ^  ♦NH  a(e))-i  (3.32) 

This  algorithm  clearly  has  similarity  to  MIM  if  (e^R  *N>  is  considered  as 
an  approximation  to  R-1  which  simply  uses  the  normalised  noise  eigenvectors 
(Gabriel  (1984)).  Alternatively.  Interpreting  (♦(g*-1  <^g)  as  an  approximation  to 
R.  the  similarity  to  PCB  becomes  apparent.  Kay  and  Demeure  (1984)  show 
that  the  MUSIC  and  PCB  estimators  are  linked  by  the  simple  relationship. 

pMU(e>  *  W  (1  -  Pp CB(0))  •  (3  33) 

for  PpcB(0)  normalised  to  a  maximum  value  of  1.  Therefore,  although  MUSIC 
Is  capable  of  producing  an  estimate  which  contains  extremely  narrow  peaks  at 
signal  directions,  whereas  PCB  results  in  a  very  flat  response,  the  'resolution' 
capability  of  the  two  algorithms  Is  Identical. 

It  should  be  noted  that,  in  the  case  of  covariance  estimates  calculated 
from  very  short  data  records,  the  noise  subspace  eigenvalues  are  not  equal, 
even  when  the  sample  is  from  a  white  Gaussian  noise  process,  MUSIC  may  be 
thought  of  as  assuming  equality  of  these  values,  perhaps  on  grounds  of 
irrelevance  or  on  the  basis  that  computer  rounding  errors  render  them 
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Inaccurate.  We  have  discovered  that  the  estimate. 

PCMU(e)  =  (a<*>H  ±NH  a(®)>_1  •  (3.34) 

often  results  In  slight  enhancement  of  the  ability  of  MUSIC  to  discriminate 
sources.  This  estimator  has  the  effect  of  increasing  the  weighting  applied  to 
those  eigenvectors  likely  to  contain  the  least  amount  of  signal.  The  EMU 
estimate  has  additional  advantages  over  MUSIC  In  that  It  more  accurately 
reflects  non  isotropic  noise  backgrounds.,  and  that  it  also  degrades  gracefully 
to  give  the  estimate  Pmlm<s>  if  no  estimate  of  signal  subspace  dimension  can 
be  formed  (Johnson  and  de  Qraaf  (1982)). 

As  demonstrated  by  Alsup  (1984)  both  MUSIC  and.,  by  implication,,  the 
PCB  filter  are  capable  of  extremely  good  discrimination  for  the  case  of 
uncorrelated  sources  when  the  noise  Is  perfectly  characterised  This  situation 
is  hardly  likely  to  arise  in  practise.  The  results  presented  here  adopt  a  more 
realistic  simulation  in  which  the  noise  Is  assumed  to  be  uncorrelated  Gaussian 
by  the  processor.,  and  the  short  data  record  is  analysed  in  an  Identity  metric. 

3.6.8.  THE  PISARENKO  MINIMUM  EIGENVECTOR  METHOD  (PME) 

The  Pisarenko  minimum  eigenvector  method  may  be  considered  as  a 
special  case  of  the  MUSIC  algorithm.,  in  which  the  noise  subspace  dimension 
equals  one.  The  eigenvector  chosen..  Is  that  associated  with  the  minimum 
eigenvalue  of  the  covariance  matrix  In  the  case  of  white  noise  corrupting  an 
Infinitely  long  data  record  (so  that  the  noise  subspace  eigenvalues  are  all 
equal) .  this  can  become  an  arbitrary  choice.  The  form  of  the  estimate  is 

PpME<e>  =  <  a(e)H  ♦mH  a(e)H  )-1  (3.35) 

Although  potentially  easier  to  compute  than  the  MUSIC  response.,  it  can  be 
seen  that.,  because  the  minimum  eigenvector  is  orthogonal  both  to  the  signal 
and  also  to  any  other  noise  related  eigenvalues.,  the  result  may  contain 
spurious  noise-generated  peaks  In  positions  unrelated  to  the  signal  directions 

3.6.9.  THE  KUMARESAN  AND  TUFTS  SIGNAL  CANCELLATION  ALGORITHMS  (KTSC) 
Many  variants  of  the  basic  signal  subspace  technique  are  discussed  in  the 

literature.  Two  methods  by  Kumaresan  and  Tufts  (1983)  for  example.,  produce 
a  single  filter  vector  which  is  effectively  a  weighted  combination  of  noise 
subspace  eigenvectors.  This  vector  is  then  used  to  formulate  the  directional 
estimates  by  comparison  with  the  array  manifold  as  before.,  but  in  some  cases 
improved  performance  over  the  MUSIC  algorithm  has  been  demonstrated 

The  two  methods  have  been  derived  by  the  authors  from  linear  predictive 
algorithms  used  for  speech  processing  (Makhoul  (1975)).  As  with  MEM  and 
many  other  algorithms.,  the  differences  are  simply  a  matter  of  notation  and 
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terminology  and  the  fundamental  approach  ot  autoregressive  modelling 
employed  is  equally  suitable  for  processing  data  of  ary  origin  which  may  be 
satisfactorily  described  by  second  moment  statistics  (van  Trees  (1968)).  given 
constram's  which  are  suitable  for  the  individual  problem  We  refer  to  the 
algorithms  as  signal  cancellation  algorithms  because  of  tho  form  of  the 
constraint  and  the  resulting  similarity  to  a  null-forming  adaptive  canceller 
equation. 

The  basic  form  of  the  equation  for  both  methods  is 

Dh  w  =  0  (3  36) 

where  the  weight  vector  is  constrained  to  be 

w  =  [  1  w2  w,  .  .  wm  1  .  (3  37) 

and  0  is  the  null  vector  In  radar  terms.,  we  wish  to  derive  an  antenna  weight 
vector  with  unity  gain  on  the  end  element,  w,  which  cancels  the  data  Since 
w,  =  1.  equation  3  36  may  be  solved  by  (along  the  first  column  of  DH.  which 
we  denote  c,  across  to  the  right-hand  side,  and  solving  for  elements  w2  to 
wm  of  w  by  forming  a  pseudo  inverse  of  the  remaining  columns  of  DH.  which 
we  denote  C.. 

w’  =  I  Wj  Wj  wm  ]  *  -1  C-1  c  (3  3S) 

Setting  the  smaller  singular  values  of  C  to  zero  ensures  that  the  solution  for  w 
is  of  minimum  length  This  is  equivalent  to  discarding  the  noise  subspace 
vectors  of  the  reduced  data  matrix,  An  excellent  and  detailed  description  of  the 
operation  of  this  algonth  n  appears  in  the  paper  by  Sibul  (1984) 

Usmg  the  singular  value  decompos  tion  of  DH.  equation  3  36  may  be 
written  in  the  form 

(Y  S  ®H)  w  =  0  (3  39) 

where  Y  and  *  are  the  appropriate  basis  vectors  Dearly  this  may  be  again 
re-written  as 

W  =  0  .  (3  40) 

Utilising  only  those  eigenvectors  which  correspond  to  the  noise  subspace.,  a 
stable  result  for  the  prediction  filter  w.  may  be  derived.,  which  can  be  shown 
to  be  of  the  form 

w  =  *n  /  (♦„  *nH  )  ..  (3  41). 

whore..  In  this  case..  <t>n  Is  the  first  row  of  the  noise  eigenvector  matrix. 
Clearly.,  this  produces  a  filter  vector  which  Is  a  complex  weighted  linear 
combination  of  the  noise  subspace  basis  vectors.,  and  compares  with  MUSIC, 
which  effectively  employs  unity  weighting  to  combine  the  eigenvectors 

Finally.,  the  spatial  estimate  is  obtained  by  correlating  the  orthogonal  filter 
vector  with  the  array  manifold.. 

pKTSC(e>  =  (  a(9)H  w  a(9)  >-1 


(3  42) 


Subsequently  we  refer  to  the  latter  form  of  the  algorithm  presented  here  as 
KTSC1.  and  the  former  as  KTSC2.  In  the  majority  of  cases  the  results 
obtained  from  these  t*.j  fo'-ms  are  Indistinguishable.  As  with  other  methods., 
the  peak  magnitudes  in  the  result  do  not  coincide  with  the  received  power  and 
so  power  levels  corresponding  to  each  of  the  directions  required  are  cerivod 
by  Interpreting  the  result  In  terms  of  point  sources  located  at  the  positions  of 
the  peaks  and  forming  the  appropriate  reduced  A  matrix.,  as  described  in 
section  3.6.3. 

The  results  produced  by  this  method  are  similar  to  those  given  by  MUSIC., 
but  with  apparently  slightly  better  noise  smoothing  and  consequently  improved 
discrimination.  However,,  because  of  the  type  of  constraint  employed  they  are 
only  ideally  suited  as  presented  here  for  regularly  spaced  and  illuminated 
arrays  of  detectors.  Care  must  be  taken  in  general  to  specify  a  suitable  form 
of  constraint  for  the  particular  problem.  An  alternative  approach,  which  wo 
have  used  with  some  success  on  non-linear  arrays  is  described  in 
section  3  8. 

3.  7.  THE  SUB-APERTURE  TECHNIQUE  (SAT) 


3.7.1.  THE  WINDOWED  DATA  MATRIX 

The  sub-aperture  technique  is  not  a  self-contained  method  of  analysis  in 
itself,  but  rather  a  useful  appendage  to  the  algorithms  discussed  above.  Its 
use  is  der.ved  from  work  on  time  series  analysis,  end  it  is  essentially  a 
method  of  forming  an  improved  covariance  matrix  estimate  from  limited  data  It 
has  important  application  to  spatial  processing  for  Its  ability  to  assist  many  of 
the  spectral  estimators  described  r.arlier  to  resolve  coherent  sources  and 
multipath,  and  also  to  convert  single  snapshot  data  into  a  suitable  format  for 
processing  by  multi-snapshot  algorithms 

The  method  consists  of  taking  the  data  vector,  and  observing  it  through  a 
small  'window'  vector  so  that  the  autocovariance  structure  of  the  series  may  be 
examined  For  every  new  observation,  we  slide  the  window  vector  further 
through  th6  data,  and  create  a  new  row  in  a  data  matrix  Thus,  if  ihe  data 
vector  is 


d  =  (dA  d,  . . ,  dn>  t3. 43) 

and  this  is  examined  with  an  r-element  window,  we  can  build  up  the  following 
matrix: 


(3  44) 
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Measuring  the  pseudo-rank  of  D  by  singular  value  decomposition  gives  the 
estimate  of  the  number  of  significant  independent  components  contributing  to 
the  data  series  if  this  is  less  than  the  minimum  dimension  of  D. 

The  forward-backward  technique  of  analysis  which  seems  popular  in  the 
field  of  speech  processing  (Makhoul  (1975) )  and  has  recently  appeared  more 
frequently  in  the  general  signal  processing  literature  (for  example  Evans  et  al 
(1982)  and  Kumaresan  and  Tufts  (1983))  creates  a  data  matrix  which,  in 
addition  to  the  'forward'  moving  window  described  above,  utilises  a  sampling  of 
the  original  series  of  data  in  the  reversed  direction.  The  reversed  data  in, 
complex  conjugated  and  included  in  the  same  matrix  as  the  forward  samples. 

It  is  clear  that  subsequent  formation  of  the  sample  covariance  matrix  involves 
extra  averaging  of  the  data  and  better  rejection  of  noise.,  particularly  in  the 
case  of  data  which  contains  sinusoidal  signal  components. 

3.  7.  g.  SINGLE  SNAPSHOTS 

If  d  is  a  data  vector  received  at  a  single  instant  in  time  from  a  uniformly 
weighted  and  spaced  linear  array  of  receiving  elements,  or  perhaps  if  it  may 
be  suitably  transformed  to  resemble  this,,  we  may  proceed  by  continuing  the 
processing  using  a  subspace  approach.  The  window  vector  is  now  the 
sub-aperture  described  by  Gabriel  (1984),  and  the  algorithm  must  be  provided 
with  an  array  manifold  for  the  sub-array  rather  than  for  the  total  number  of 
elements.  This  is  simply  computed  as  a  sub-section  of  the  original 
translationally  invariant  manifold.  Having  done  this,  a  startling  improvement  in 
the  ability  to  discriminate  closely  spaced  point  emitters  is  possible. 

3.  7.  3.  MULTIPLE  SNAPSHOTS  :  COHERENT  SOURCES 

Sub-aperture  processing  may  also  be  beneficially  applied  to  multiple 
snapshot  data.,  particularly  for  Its  effect  on  data  from  coherent  sources.  In  this 
case  It  Is  perhaps  easier  to  think  of  using  an  (r  x  m)  rectangular  sampling 
matrix  K  of  the  form., 

IJ.  =  [  0  I  I  I  0  1  ..  (3.  45) 

where  the  subscript  I  denotes  the  column  in  which  a  1  first  appears,  to  extract 
a  principal  submatrix..  Rj..  of  the  original  sample  covariance  matrix., 

Rj  =  jC  R  (3.48) 

Averaging  the  submatrices  produced  in  this  manner  then  has  rho  red  resuh 
(Evans  et  al  (1982) ) . 

For  the  case  of  data  originating  from  coherent  sources.,  the  improvement 
attainable  with  this  approach  may  be  understood  as  deriving  from  the  'internal 
averaging'  of  the  normally  computed  correlation  function  (Evans  et  al  (1982)) 
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which  causes  the  non-stationary  terms  to  tend  toward  zero.  Alternatively. 
Gabriel  (1982)  suggests  that  the  improvement  may  be  interpreted  as  arising 
from  the  movement  of  the  sub-array  phase  center  with  respect  to  the  targets, 
which  introduces  an  effective  doppler  phase  shift  between  the  sources. 

3.7.4.  TIME  SERIES  PROCESSING 

Clearly,  such  processing  can  have  wide  applicability.,  and  may  perhaps  be 
used  on  data  from  domains  other  than  the  spatial  domain.  It  was  suggested  in 
section  3.6.3.  for  example.,  that  more  information  about  the  nature  of  the 
received  signals  could  be  deduced  by  investigating  the  time  series  component 
of  the  data  in  isolation.  This  is  a  potentially  useful  post-processing  stage  after 
reducing  the  number  of  time  series  to  be  examined  by  means  of  the  spatial 
spectral  estimators  described  previously  Taking  an  individual  time  series  from 
the  matrix  T  (equation  3.20).  corresponding  to  a  particular  direction  identified 
as  worthy  of  further  investigation,  it  is  possible  to  estimate  the  number  of 
decorrelated  signals  present  at  that  direction  by  forming  a  covariance  estimate 
using  the  sub-aperture  technique  and  processing  the  result  by  the  same 
methods  used  for  spatial  data  (Zlegenbein  (1979)).  In  fact,  if  a  library  of 
expected  time  domain  signals  is  available,  the  individual  waveforms  which  are 
present  may  be  estimated.  In  the  absence  of  such  a  manifold.,  the  eigenvalues 
of  the  covariance  matrix.  (D  DH) .  may  simply  be  used  to  give  an  estimate  of 
the  relative  powers  of  the  different  component  modes,  given  sufficient  samples 
(van  Trees  (1968))  The  corresponding  problem  of  estimating  source 
locations.,  given  the  emission  frequencies  may  also  be  tackled  in  this  way 

The  benefits  of  using  such  a  sequential  approach  to  the  processing  are 
illustrated  by  many  of  the  examples  presented  in  Section  4.  If  data  has  been 
gathered  in  more  than  two  domains  (for  example  along  two  spatial  axes  and  in 
time)  it  is  clear  that  this  sequential  strategy  may  be  extended  to  process  It  in 
a  simple  manner. 

3.7.5.  SUMMARY 

It  can  be  seen  that  the  sub-aperture  technique  enables  the  detection  of  a 
larger  numoer  of  independent  sources  than  the  rank  of  the  usual  data  matrix,, 
limited  as  It  is  by  the  minimum  dimension  of  the  matrix,  might  initially 
suggest.  We  have  shown,  for  example.,  that  it  is  possible  both  to  locate 
coherent  emitters  which  are  sufficiently  spatially  uncoi  related  and  also  the 
frequencies  of  spatially  correlated  targets  It  they  are  decorrelated  in  the  time 
domain.  This  suggests  that  the  Introduction  of  data  from  additional 
measurement  domains  In  which  targets  are  decorrelated..  will  not  only  allow 
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more  detailed  characterisations  of  sources  already  detected.,  but  will  also 
Increase  the  total  number  of  sources  which  the  processor  is  able  to  detect.  If 
a  sequential  processing  strategy  is  adopted,  tackling  the  problem  which  has 
the  smallest  dimension  first  enables  the  remaining  load  to  be  accordingly 
reduced.  However.,  if  targets  are  not  easiiy  discriminated  in  that  domain,  some 
information  may  be  lost  or  inaccurately  recovered,  as  will  be  discussed  in 
section  3  9  2,  and  demonstrated  in  section  4.2.8. 

3.8.  REDUCING  THE  COMPUTATIONAL  LOAD 

3.8.1.  COMPARISON  OF  THE  ALGORITHMS 

If  computational  simplicity  is  the  major  criterion  in  choosing  a  suitable 
algorithm  for  signal  processing,  then  clearly  the  optimum  choice  is  either  the 
BSA  or  the  PIM  algorithm.  Both  Involve  pre-computation  of  the  filter  matrix, 
and  so  data  processing  can  take  place  on  line,  acting  on  each  snapshot  as  a 

simple  matrix-vector  multiplication,  or  on  blocks  of  data  in  the  form  of  the 

covariance  matrix  estimate.  If  the  enhanced  discrimination  of  the  other 
algorithms  presented  is  desired.  however,  the  computational  load  is 
significantly  increased  and  calculations  may  have  to  be  carried  out  off  line 

The  true  maximum  likelihood  method,  although  optimal  from  the  point  of 
view  of  estimation  theory  (van  Trees  (1968)).  rapidly  becomes  computationally 
intractable  as  the  number  of  assumed  wavefronts  is  increased,  since  It  involves 
both  matrix  Inversions  and  a  simultaneous  multi-parameter  search  for  each  new 
sample  covariance  matrix.  It  is  more  common  to  find  it  used,  as  we  have 

done  for  the  current  work,  in  the  form  given  in  section  3.6  4  which  has  been 

derived  on  the  assumption  of  a  single  plane  wave  arriving  at  the  antenna. 

If  a  full  singular  value  analysis  of  the  data  matrix  is  needed  for  the  simpler 
KTSC  algorithm,  for  example,  this  might  be  accomplished  by  diagonalisatlon  of 
the  covariance  matrix  using  standard  routines.  Householder  tridiagonalisation  of 
the  matrix,  followed  by  calculation  of  all  the  eigenvalues  and  eigenvectors 
using  a  QR  routine  (Wilkinson  (1965),  Wilkinson  and  Reinsch  (1971))  would 
involve  2ms/3  multiplications  and  m  square  roots  for  the  first  stage,  and  a 
further  4m*  multiplications  with  m-1  square  roots  for  each  iteration  of  the  QR 
algorithm  (generally  the  number  of  iterations  required  is  of  the  order  of  two 
per  eigenvalue  for  the  tridiagonalised  matrix).-  in  addition,  the  formation  of 
rank-reduced  pseudo-inverses  and  matrix  approximations  requires  yet  more 
processing  for  each  new  block  of  data.  If  It  Is  to  be  carried  out  on  line. 

Because  of  the  usual  ms  dependency  of  the  processing  load,  the 
requirement  escalates  rapidly  with  the  size  of  the  problem,  and  so  from  this 
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point  of  view  algorithms  which  reduce  the  size  of  the  matrix  to  be  dealt  with 
are  potentially  very  valuable.  An  example  of  such  an  algorithm  is  the 
sub-aperture  technique,  which  reduces  the  covariance  matrix  to  the  dimension 
of  the  window  used.,  at  the  expense  however  of  a  possibly  arbitrary  reduction 
In  the  number  of  available  degrees  of  freedom  for  modelling  the  targets. 

It  is  worth  noting,  however,  that  for  problems  involving  only  a  small 
number  of  sensors.,  the  eigen  analysis  may  actually  require  less  computation 
than  is  then  needed  for  the  location  of  the  sources.,  if  this  involves  searching 
a  large  number  of  points  Instead  of  exhaustive  searching,  more  efficient 
algorithms  must  be  employed  For  searches  Involving  more  than  one  data 
domain,  our  sequential  processing  approach.,  which  uses  the  results  of  a 
search  in  each  domain  as  support  for  a  reduced  search  in  the  next  domain, 
has  proved  to  be  effective  in  reducing  the  complexity  of  the  computation 
(Clarke  et  al  (1985)) . 

3.  8.  2.  THE  COEFFICIENTS  METHOD 

We  know  from  section  3.6.2  that  the  system  response  matrix..  A,  may  be 
a  less  than  full  rank  transform  and  that  our  data..  D,  will  therefore  generally 
contain  redundant  or  noisy  components.  Analysis  of  the  data  in  terms  of  the 
aporopriate  basis  vectors  of  A.  and  rejection  of  components  at  or  below  the 
noise  floor  allows  us  to  represent  the  useful  data  in  a  more  compact  form 

As  a  result  of  the  singular  transformation.. 

A  F  +  N  =  D  .  (3  47) 

D  is  a  matrix  of  dimensions  (m  x  NpK.  where  m  is  the  number  of  elements  in 
the  array,  and  Np  Is  the  number  of  snapshots  In  the  block  of  data  to  be 
processed . 

Premultiplying  both  sides  of  this  equation  by  the  transposed  set  of 
left-hand  singular  vectors  of  A  gives 

Uh(AF  +  N)  =  UH  D  =  C  (3  48) 

In  general,  the  size  of  the  matrix  U  may  be  reduced  by  examining  the 
singular  values  of  A  in  relation  to  the  estimated  noise  floor  This  number  will 
be  less  than  m,  and  may  be  denoted  k  Thus,  the  data  is  now  represented  by 
Its  projection  onto  the  “significant"  basis  vectors  of  the  array  manifold,  and  is 
stored  in  a  matrix.  C.  of  dimensions  (k  x  Np).  This  has  allowed  a  reduction 
in  the  size  of  the  effective  data  array  to  a  minimum  size  consistent  with 
retaining  all  potentially  useful  degrees  of  freedom. 

It  is  now  possible  to  use  the  matrix,  C.  as  the  input  to  a  signal-subspace 
processor.  However,  the  assumed  array  manifold  must  also  be  projected  onto 
the  same  basis  vectors,  and  becomes  (U*'*  A).  The  vector  or  set  of  vectors 
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orthogonal  to  the  signal  subspace  which  the  algorithms  produce  can  be 
compared  against  the  columns  cf  this  transformed  manifold  in  order  to  identify 
source  directions,  in  this  way.,  it  is  possible  to  significantly  reduce  the 
computational  load  in  poor  signal  to  noise  environments,  when  the  rank  of  A  is 
effectively  very  much  lower  than  the  number  of  detector  elements. 

We  refer  to  this  procedure  as  the  coefficients  method  since  the  new  data 
matrix  may  be  thought  of  as  consisting  of  coefficients  of  the  known  array 
manifold  basis  vectors. 

The  projection  of  data  onto  an  alternative  basis  set  of  vectors  in  this  way. 
which  may  be  interpieted  as  an  additional  linear  transformation  of  the  data, 
has  proved  of  some  use  in  processing  the  data  from  nonlinear  arrays  of 
receivers.  In  the  ca^a  ct  a  circular  array  of  detectors,  operating  only  in  the 
plane  of  the  array  and  with  the  rear  semi-circle  assumed  to  be  shielded  from 
a  source  at  a  particular  angle.,  the  data  matrix  contains  a  large  number  of 
zero  entries.  If  the  KTSC  as  described  is  understood  as  producing  a  sidelobe 
canceller  type  of  behaviour.,  then  it  can  be  seen  that  this  format  of  data  matrix 
leads  to  an  insoluble  equation  if  we  cancel  against  a  master  element  which  is 
hidden  from  the  main  emitters.  Transforming  the  data  by  the  mapping  just 
described  removes  such  difficulties  in  addition  to  decreasing  the  time  spent  in 
computation. 

3.9.  ITERATIVE  PROCESSING 

3.9.1.  ITERATING  THE  PSEUDO-INVERSE 

The  pseudo-inverse  method  described  in  section  3.6.2  depends  partly  on 

the  definition  of  limits  to  the  angular  extent  of  sources  for  its  ability  to 

enhance  the  discrimination  of  close  emitters.;  The  proportionate  reduction  in 
the  width  of  the  system  impulse  response  over  that  of  the  BSA  matched  filter 
increses  slightly  as  the  defined  angle  of  view  Is  restricted  for  a  given  signal  to 
noise  level.  Thus  it  is  possible  to  use  the  result  of  a  reconstruct-on  which 
assumed  a  wide  angle  of  view  to  provide  a  form  of  prior  information  for  an 

iterative  form  of  processing  of  the  data  In  the  hope  that  this  will  lead  to 

improved  discrimination.. 

Such  processing  might  be  carried  out  by  equating  values  in  the  output  of 
the  initial  estimate  which  fall  below  a  threshold  level,  e.  to  that  same  value., 
and  all  remaining  values  to  unity.  The  resulting  function  would  then  be  used  to 
weight  the  corresponding  columns  of  the  A  matrix  before  forming  a  new 
inverse.  This  will  require  intensive  computation  for  the  processing  of  each 
block  of  data,  thus  removing  one  of  the  major  advantages  of  PIM.  It  will  be 
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shown  by  results  presented  in  section  4  that,  although  this  method  ensures 
that  sources  resulting  In  measured  amplitudes  which  fall  above  e  are  not 
further  attenuated.,  any  low  amplitude  sources,,  which  might  be  easily  detected 
by  alte  native  algorithms.,  could  easily  fall  below  a  simple  threshold  based  on 
the  sideiobe  levels  of  PIM  and  thus  would  be  strongly  suppressed  The  ensuing 
reconstructions  would  be  seriously  in  error.  This  approach  also  results  in 
spectra  which  contain  sharp  discontinuities  at  the  switching  boundaries. 

An  alternative  approach  could  perhaps  involve  setting  the  below  threshold 
regions  of  the  reconstruction  to  e,  and  all  otner  values  to  a  level  reflecting  the 
magnitude  of  the  latest  reconstruction  at  that  point.  It  is  not  clear  from  our 
results  that  this  gives  any  significant  gain  in  discrimination  although  it  may 
sometimes  be  used  to  beautify  the  result.  As  we  have  found  m  practice,  this 
method  also  tends  to  suppress  low  amplitude  sources  still  furthor  with  each 
successive  iteration  A  simitar  result  would  obtain  from  applying  such  a 
weighting  to  an  iterated  BSA  and  would  be  considerably  easier  to  evaluate 

If  we  interpret  the  original  A  matri>  as  Implying  a  'top  hat'  weighting  on 
the  space  of  all  possible  solutions  He  unity  weighting  to  di- actions  within  the 
'top  hat'  and  zero  weighting  to  all  other  directions) .  the  methods  described  in 
this  section  may  be  interpreted  as  placing  a  further  emphasis  on  certain 
directions  within  the  space  of  solutions  so  that  the  weighting  function  reflects 
the  shape  of  the  reconstruction  at  each  iteration. 

These  approaches  seem  to  have  much  In  common  with  the  yet  more  easily 
computed  practice  of  squaring  the  magnitude  of  results  in  order  to  sharpen 
already  noticeable  peaks  and  to  suppress  sideiobe  effects  This  too  gives 
additional  weighting  to  particular  directions  in  the  space  of  possible  solutions., 
as  Is  borne  out  by  results  which  we  present  in  Section  4  2  9 

3.9.2.  USING  SUBSPACE  ALGORITHMS  AS  PRIOR  KNOWLEDGE 

Since  a  computer  cannot  store  the  infinite  number  of  vectors  required  for 
a  complete  description  of  a  continuous  space.,  the  array  manifold.  A  is 
generally  stored  in  discrete  form.  For  this  reason  the  steering  vectors  used  in 
the  initial  estimates  may  be  slightly  mismatched  to  the  actual  source 
directions  It  is.,  however.,  a  simple  matter  to  improve  the  estimates  by 
iterative  interpolation  of  steering  vectors  to  search  for  maximum  peak  amplitude 
or  null  depth  in  the  regions  identified  by  the  initial  coarse  processing 

Having  obtained  estimates  of  the  number  of  signals  and  their  directions  of 
arrival  by  means  of  the  signal  subspace  algorithms,  these  can  be  interpreted 
as  prior  knowledge  for  further  processing  designed  to  extract  greater  knowledge 
of  the  target  characteristics.,  as  described  In  sections  3  6  3  and  3  7  The 
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precise  directional  estimates  obtained  from  signal  subspace  algorithms  may 
also  be  taken  as  prior  knowledge  for  the  choice  of  additional  steering  vectors 
to  be  used  in  the  array  manifold..  A.  as  used  for  the  PIM  processing  described 
in  section  3.6  2 

It  is  also  worth  noting  here  that  the  order  in  which  the  data  domains  are 
investigated  by  a  sequential  implementation  of  the  signal  subspace  algorithms 
is  of  particular  importance  in  certain  cases  Singular  value  analysis  of  3  data 
matrix  may.,  for  example.,  suggest  the  presence  of  two  incoherent  sources 
wNch  are  not  subsequently  discriminated  by  orthogonal  filtering  in  one  of  the 
domains.  For  example.,  two  sources  which  are  highly  correlated  in  the  spatial 
frequency  domain,  but  which  are  decorrelated  in  time  may  not  be  discriminated 
spatially.,  despite  the  data  matrix  clearly  having  a  pseudo  rank  of  two.  in  this 
case.,  the  sources  ma>  a  discriminated  first  of  all  by  examining  the  time 
domain  singular  vectors.  The  resulting  frequency  information  can  then  be  used 
to  extract  corresponding  independent  data  series  from  the  data  matrix, 
enabling  the  emitter  directions  to  be  derived  with  potentially  greater  accuracy 
and  ease. 

Finally.,  it  should  be  noted  that,  whatever  the  validity  of  interpretation  of 
the  data  in  terms  of  point  targets  (this  is  unlikely  to  be  a  suitable  model  for 
all  possible  radar  targets  for  example)  .,  it  is  undeniably  a  very  much  more 
compact  representation  than  the  data  matrix  itself.  The  exercise  may  therefore 
be  viewed  in  its  alternative  role  as  simply  resulting  in  data  compaction. 
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4.  COMPUTER  SIMULATIONS 


4.1.  THE  SIMULATION  FACILITY 


A  computer  simulation  facility  has  been  developed  ir.  order  to  test  and 
compare  the  performance  of  the  techniques  described.-  it  enables  the 
characterisation  of  co-linear  and  polygonal  aperture  plane  and  cc-llnear  focal 
plane  arrays  of  detectors  for  the  cases  of  both  Ideal  and  perturbed  weight 
vectors.  Array  manifolds  may  be  analysed  In  terms  of  their  singular  value 
spectra,  and  the  relevant  data  stored  for  future  use  on  a  library  disc  A  menu 
driven  structure  allows  the  creation  of  target  scenarios  either  for  immediate  use 
or  for  storage  on  the  library  disc  and  experiments  have  been  carried  out  using 
simulations  of  both  regularly  and  randomly  sampled  sine  wave  emitters.,  and 
also  of  narrow  band  f.  m.  sources  Data  generated  by  running  such 
simulations  with  the  addition  of  spatially  uncorrelated  Gaussian  noise  may  be 
processed  by  any  of  the  algorit.ims  described  in  Sect-ons  3. 6  and  3  7  of  this 
report  and  the  results  may  be  d*rectly  compared  by  running  each  technique 
with  identical  data  samples  A  rsoresentative  selection  of  the  results  obtained 
Is  presented. 

4.g.  DISCUSSION  OF  RESULTS 

4.  2.  1 ,  NOTATION 

Our  results  are  generally  presented  in  the  torm  of  plots  of  power  against 
number  of  cycles.  The  number  of  cycles  relates  to  the  spatial  or  temporal 
frequency  of  interest.,  as  discussed  in  Section  3.3.2.  Vertical  lines  surmounted 
by  diamond  shapes  indicate  the  signal  frequencies  and  powers  used  to  create 
the  data  simulations.  Vertical  lines  surmounted  by  norizomal  bars  indicate  the 
signal  frequencies  and  total  powers  as  perceived  by  the  signal  subspace 
algorithms,  taking  peak  centres  as  the  measure  of  frequency 

4.2.2.  SENSITIVITY  OF  PIM  AND  BSA  TO  PRIOR  ASSUMPTIONS 

The  capability  of  the  PIM  beamformer  to  discriminate  the  presence  of  even 
a  single  source  of  power  depends  on  the  formation  of  a  suitably  supported 
matrix  inverse  The  angular  limits  to  the  field  of  view  must  be  defined  befc'e 
singular  value  decomposition  can  be  earned  out.  and  knowledge  ot  the  signal 
to  noise  ratio  Is  required  for  the  formation  ot  a  stable  matrix  Inverse,  as  will 
be  shown.  Strictly,  the  signal  to  noise  ratio  at  each  component  singular 
frequency  must  be  known  in  order  to  find  the  best  possible  inverse,  but  our 
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FiO  4.  la.  Raaponaa  or  PIM  ana  QSA  to  a  alngia  amlttar.  Signal  to  nolaa  ratio 
la  130  dB  Tha  array  la  a  20  alamant  linaar  array,  and  20  anapahots  of  data 
hava  baan  procaaaad. 


t 

Fig  a  10  Raaponaa  of  PIM  and  BSA  to  a  alngia  amlttar.  Signal  to  nolaa  la 
SO  dB.  Out  tha  PIM  raconatructlon  haa  aaaumad  a  noiaa  (aval  of  only  -130  dB. 
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Fig  4  It.  Response  o!  PlM  and  BSA  to  a  single  emitter  Signal  to  nolsa  is 
0  dB.  Tha  beamwiptne  of  tha  two  methods  ara  virtually  Identical 
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Fig  4  2a  Response  of  PlM  and  BSA  to  two  equal  power  sources,  separated  by 
0  8  conventional  beamwidths  PlM  has  achieved  'classical''  resolution. 


simulations  have  assumed  this  information  not  to  be  available. 

Pig  4.1  shows  the  response  of  the  PIM  and  6SA  beamformers  to  a  single 
emitter.,  located  broadside  to  a  co-linear  array  of  Shannon  number  three.,  for 
a  wide  range  of  signal  to  noise  ratios.  We  assume  initially  (Fig  4.  la)  that  the 
singular  value  spectrum  of  the  array  manifold  has  been  truncated  at  a  level 
corresponding  to  the  presence  of  a  noise  floor  at  the  unrealistically  low  level 
of  -130  dB.  The  response  In  this  virtually  noise  free  environment  displays  the 
narrower  beamshape  of  PIM  over  that  of  BSA.  Clearly,,  even  at  this  S/N  the 
difference  in  width  is  not  great. 

Fig  4.1b  shows  the  result  of  an  incorrect  assumption  about  the  level  of 
additive  noise  in  the  formation  of  the  pseudo  Inverse.  The  plot  of  PpiM(e)  has 
been  computed  assuming  a  noise  level  of  -130  dB,  but  with  a  simulation 
incorporating  noise  at  -80  dB,  The  result  is  a  meaningless  amplification  of  the 
noise  components  of  the  data.  Fig  4,1c  shows  the  result  obtained  by  using  a 
more  suitable  truncation  of  the  singular  values.  It  is  perhaps  worth  noting  thai. 
even  at  this  still  very  high  signal  to  noise  ratio,  we  have  only  needed  ten  of 
the  original  twenty  array  manifold  singular  values  for  the  reconstruction.  At  the 
very  much  more  realistic  S/N  of  10  dB  (Fig  4,  Id)  and  using  only  the  six  most 
significant  singular  values,  the  beamwidth  of  the  PIM  estimator  is  appreciably 
broader  than  before,  and  for  a  S/N  of  0  dB,  using  three  singular  values, 
becomes  indistinguishable  from  that  of  the  BSA  (Fig  4.  1e). 

Since  the  enhancements  obtainable  from  PIM  are  only  clearly  observed  in 
conditions  of  high  signal  to  noise,  ensuing  simulations  in  this  section  will  be 
demonstrated  in  relatively  noise  free  environments. 

First  of  ail.  Fig  4.2a  demonstrates  the  classical  'resolution'  problem,  w.th 
two  emitters  cf  50  dB,  separated  by  0.8  cycles  across  the  aperture  In  the 
conventional  sense,  tahing  only  the  major  peaks  as  corresponding  to  resolved 
emitter  positions,  PIM  demonstrates  superior  resolution  capability  when 
compared  to  the  BSA.  However.  Fig  4.2b  demonstrates  very  clearly  the 
difficulty  with  this  conventional  view  of  resolution,  when  the  second  emitter  is 
of  lower  power  than  the  first  The  only  indication  of  the  presence  of  the 
second  target,  as  can  also  be  seen  in  the  BSA  plot  of  Fig  4.2a,  is  in  the 
asymmetry  of  the  result  about  the  main  lobe  of  the  response.  We  know  this 
because  we  have  prior  knowledge  about  the  expected  shape  of  the  response  to 
a  single  emitter  that  we  expect  a  sin(x)/x  response  from  the  BSA,  for 
example.  Given  at  least  two  independent  snapshots  of  the  data,  a  signal 
subspace  algorithm  such  as  MUSIC  would  hanc  this  example  with  ease, 
producing  sharp  clearly  defined  peaks  corresponding  to  the  source  directions 

Moving  the  lower  power  interfering  source  outside  the  limiting  angle  of  view 
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defined  *o r  P!M,.  we  obtain  Fig  4.3a  Whilst  the  BSA  beamformer  remains 
unperturusd  by  this,  the  PIM  plot  demonstrates  once  again  the  extreme 
sensitivity  of  this  algorithm  to  the  validity  of  the  prior  information  required  for 
the  formation  of  the  matrix  pseudo  inverse.  The  directional  estimate  produced 
by  MUSIC  for  the  central  emitter  (Fig  4.3b)  is  extremely  precise.,  but  the 
corresponding  estimate  of  signal  power.,  since  it  utilises  a  pseudo  inverse 
supported  with  only  a  single  delta  function,  is  in  error.,  as  demonstrated  by 
the  high  power  attributed  to  a  low  level  spurious  peak. 

We  now  demonstrate  how  KTSC.  or  a  similar  method  such  as  MUSIC.,  may 
be  used  to  gain  prior  knowledge  lor  the  support  of  the  PIM  pseudo  inverse 
required  by  the  above  problem.  Since  the  weight  vector  derived  by  the  KTSC 
algorithm  is  derived  from  the  data  alone,  the  assumption  of  a  narrow  field  of 
view  required  by  PIM  is  not  neccessary.  and  as  we  saw  in  Fig  4  3b  the 
accuracy  of  the  directional  estimate  is  therefore  not  susceptible  to  such 
assumptions  Simply  by  scanning  a  much  wider  angle  of  view  we  may  obtain 
the  spatial  frequency  of  the  Interfering  source.  This  frequency  corresponds  to 
e  column  of  a  much  broader  A  matrix  than  was  used  for  the  PIM  estima'e.  If 
we  simply  add  this  single  column  to  the  original  A  matrix,  ana  re-compute  the 
pseudo  inverse,,  the  normalised  result  shown  in  Fig  4.3e  is  obtained.  The 
deleterious  effect  of  the  interfering  source  has  been  all  but  completely 
removed  For  as  long  as  this  source  remains  spatially  stationary.,  the  new 

pseudo  inverse  may  be  used  to  provide  a  rapid  spectral  estimate  by  simple 
matrix  multiplication  on  each  block  of  data  received. 

4.2.3.  UTILISATION  OF  THE  TIME  DOMAIN  INFORMATION 

If  the  sources  we  are  attempting  to  detect  emit  sinusoidal  waves  of  slightly 
different  frequencies,  it  Is  Clear  that  the  spatial  power  distribution  simply 
consists  of  a  three  dimensional  interference  pattern  which  shifts  as  the  sources 
appear,  on  a  snapshot  by  snapshot  basis.,  to  change  their  relative  phases 
This  is  perceived  across  the  face  of  the  linear  array  of  antennae  by  a 
sinusoidal  variation  of  signal  amplitude  sweeping  across  the  array  As  this 

occurs.,  if  we  perform  a  PIM  type  reconstruction  of  each  snapshot  as  it  is 

received.,  we  see  that  the  result  may  be  dependent  on  that  relative  phase 

between  the  two  targets.  V.'ltcr  the  pair  appears  in  antiphase  as  observeo  at 
the  antenna  phase  centre.,  this  is  manifested  as  a  null  in  the  power  distribution 
across  the  elements.,  and  a  maximised  power  gradient  across  the  arr?  /  This 
enables  the  beamformer  to  more  easily  discriminate  the  presence  of  two 
targets.  Conversely,  when  the  sources  appear  to  be  in  phase.,  the  power  at 
the  army  phase  centre  :s  a  maximum.,  and  the  power  gradient  in  the  local 
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Fig  4.4«.  A  snapshot  by  snapshot  reconstruction  using  PtW  to  demenstrat#  the 
effect  of  relative  largat  phase  for  two  aiowly  beating  aourcam. 


Fig  4.4b  The  result  of  PIM  and  BSA  operating  on  the  covariance  matrix 
formed  from  the  sixty  snapshots  cf  Fig  4  4a  Dhta  fa  taken  from  a  twenty 
element  Una  array  of  sensors  spaced  by  x/2.  The  tower  curve  la  the  function 
1/t  t-Ppjj*(e> )  to  show  the  positions  of  the  major  peaks  of 


space  Is  a  minimum.  Such  a  situation  provides  a  difficult  problem  for  the 
beamformer.  Fig  4.4a  shows  a  series  of  such  single  snapshot  PiM 
reconstructions,  with  power  plotted  along  a  vertical  linear  axis.  Such  a  plot 
tells  us  that  there  must  be  more  than  one  point  emitter  in  the  object  space, 
but  not  where  the  sources  might  be  located,  or  the  possible  source  powers. 
We  will  now  attempt  to  derive  this  information  from  our  multiple  snapshot 
algorithms. 

Combining  the  snapshots  through  the  formation  of  a  sample  covariance 
matrix,  we  can  plot  the  multiple  snapshot  results  from  the  PIM  and  BSA 
processors  for  this  example  (Fig  4.4b).  Also  plotted  In  this  figure  is  the 
function  1/<l-f(x>).  where  f(x)  is  the  normalised  PIM  result.  This  shows  the 
positions  of  the  main  PIM  peaks  which  have  just  beon  resolved,  and  indicates 
that  ihey  do  not  correspond  to  the  directions  from  which  the  signals 
originated.  Fig  4.4c  demonstrates  that  this  case  is  handled  with  ease  by 
MUSIC,  which  has  determined  both  signal  powers  and  directions  with  great 
accuracy.  We  conclude  from  this  that  examination  of  the  data  matrix  rather 
than  the  array  manifold  enables  superior  resolution  of  the  type  of  problem 
Involving  docorrelated  point  emitters,  through  better  utilisation  of  the  additional 
Information  available  from  the  time  domain  decorrelation  (a  single  point  source 
would  not  result  In  the  paltorn  shown  in  Fig  4.4a). 

4.2.4.  COMPARISON  OF  ALL  ALGORITHMS  ON  SAME  DATA 

We  have  carried  out  numerous  comparative  simulations  of  the  algorithms 
coverod  by  this  report,  and  present  here,  as  an  examplo.  P(e)  plots  for  each 
of  the  processors  acting  on  data  from  a  simple  scenario  involving  three 
emitters  (clg  4.5).  The  simulation  is  of  a  twenty  element  half  wavelangth 
spaced  collnear  array,  with  a  Shannon  number  of  six.  Data  has  boen  sampled 
at  a  very  low  rate  to  simulate  the  undersampling  likely  to  occur  in  a  real 
digital  to  analogue  conversion  system,  and  has  thus  allowed  the  sources  ample 
time  to  decorrelate.  The  sources  consist  of  two  12  dB  sinusoidal  emitters  In 
close  proximity  to  a  60  dB  random  phase  jamming  signal.  In  all  but  the  BSA 
and  PiM  cases  we  have  derived  Ume  series  from  the  data  which  correspond  to 
each  of  the  peaks  Identified,  and  analysed  them  with  a  ten  element  moving 
sub-aperture,  as  described  In  section  3.7.4.  The  eigenvalues  of  these 
analyses  have  been  taken  as  estimates  of  the  power  of  independent  sources, 
end  plotted  as  arrow  heads  on  the  lines  corresponding  to  the  respective 
directions.  This  enables  us  to  examine  the  amount  of  leakage  of  the  jamming 
signal  into  the  Identified  signal  directions. 

Referring  to  Fig  4.5,  plot  a  shows  the  response  of  the  MUSIC  processor 
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Mg  4  5t  The  data  of  Fig  4  5a  procassad  by  TNA 
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Fig  4  5g  Tna  results  obtained  from  processing  the  same  data  as  the  previous 
figure,  using  the  PCB  algorithm  (shown  by  the  upper  curve)  The  lower  curve 
*how*  the  function  l/(  1~PpcB<*>  >  to  demonstrate  Its  similarity  to  the  MUSIC 
curve  of  Fig  4.5a 
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to  this  simulation.  The  two  signal  directions  are  accurately  determined.; 
although  only  lust.,  as  Indicated  by  the  shallowness  of  the  dip  between  the 
relevant  peaks  The  area  away  from  the  signal  sources  is  depicted  as  being 
very  flat,  as  a  result  of  the  averaging  of  the  seventeen  equally  weighted  noise 
eigenvectors  used  by  MUSIC  In  this  example.  The  three  very  shallow  peaks 
Identified  in  this  region  are  found  to  be  associated  with  total  power  levels 
significantly  below  the  noise  level.,  and  can  safely  be  ignored.  The  cluster  of 
high  eigenvalues  derived  from  the  time  series  in  the  direction  of  the  jammer 
Indicates  the  presence  of  a  locally  very  high  noise  floor.,  in  each  of  the  signal 
directions.,  a  single  eigenvalue  stands  out  higher  than  the  rest,  indicating  the 
presence  of  the  sinusoidal  emitters  in  a  heavily  rejected  noise  background  from 
the  jammer. 

Plot  b  depicts  the  result  of  processing  the  same  data  by  the  KTSC1 
algorithm.  Utilising  exactly  the  same  eigenvectors  as  did  MUSIC.,  we  obtain  a 
different  result.  The  signal  sources  have  resulted  in  two  clearer  peaks., 
although  the  particular  sample  of  noise  present  in  the  data  has  resulted  in 
slightly  poorer  estimates  of  direction  in  this  case.  The  background  ripple  is 
more  pronounced  than  in  MUSIC.,  although  such  peaks  may  still  be  rejected  on 
the  basis  of  the  low  power  levels  detected  there.  Interpreting  the  KTSC1  weight 
vector  as  a  polynomial  in  the  z-plane  (Kumaresan  and  Tufts  (1983)).  plot  c 
shows  the  stable  regular  distribution  of  zeroes  around  the  interior  of  the  unit 
circle  The  signal  zeroes  lie  almost  on  the  circle. 

Plot  d  shows  the  result  obtained  using  MEM  on  this  data  sample.  Although 
both  signal  and  |ammer  directions  have  been  accurately  determined.,  note  the 
typical  presence  of  pronounced  peaks  at  other  directions.  The  proximity  of  one 
of  these  peaks  to  the  signal  frequencies,,  and  the  consequent  relatively  high 
power  level  associated  with  it  has  resulted  in  higher  levels  of  noise  leakage 
from  the  jammer  into  the  signal  directions  as  perceived  through  the  incorrectly 
supported  pseudo  inverse.  A  simple  iteration  which  removed  the  spurious  peaks 
from  the  support  for  this  operation  would  result  in  improved  power  estimation 

MLM  (plot  e)  has  tended  to  average  out  the  major  spurious  peak,  at  the 
expense  n  this  case  of  less  accurate  determination  of  the  signal  directions. 
The  result  of  TNA  processing,  shown  in  plot  f  has.  not  surprisingly., 
similarities  to  both  MEM  and  MLM. 

Plot  g  shows  both  the  PCB  estimate  and  the  function  l/(l-f(x>>  derived 
from  equating  f(x>  to  the  normalised  matched  filter  estimate.  This  result  clearly 
emphasises  the  fact  that  the  visual  appearance  of  an  estimate  is  not  a  reliable 
indicator  of  its  ability  to  resolve  targets  in  the  sense  defined  by  Rayleigh.  The 
inter-peak  dips  of  the  matched  filter  form  of  this  estimate  are  of  only  a 


fraction  of  a  dB.  and  yet  they  are  stable  because  of  the  removal  of  noise  in 
the  rank  reduced  sample  covariance  estimate.  The  major  error  evident  in  this 
plot  (the  hopelessly  Inaccurate  estimate  of  signal  powers)  has  arisen  not  as  a 
result  of  a  deficiency  of  the  PCB  algorithm,  but  through  the  inability  of  the 
simple  peak  locating  routine  employed  to  accurately  find  the  centre  o:  the 
extremely  broad  response.  This  has  caused  a  slight  inaccuracy  in  the  pseudo 

inverse  support  for  the  important  iammer  direction,  and  hence  significant 

leakage  of  this  relatively  high  power  signal  into  the  signal  time  series.  This 

only  serves  to  demonstrate.,  however.,  the  high  degree  of  cancellation  which 
may  possibly  be  achieved  by  this  form  of  processing  if  properly  executed. 

PEMU(e).  shown  In  plot  h  Is  very  similar  to  the  MUSIC  result  in  this 
case.,  but  shows  by  the  presence  of  more  prcnounced  peaks  in  the  background 
floor  that  the  flat  floor  of  MUSIC  is  partially  caused  by  the  assignment  of 
uniform  weighting  to  each  of  the  eigenvectors  in  the  noise  subspace.  The  plot 
obtained  .'rom  the  PME  estimator  (plot  I)  demonstrates  that  the  flat  floor  of 
MUSIC  is  silso  the  result  of  averaging  many  noise  eigenvectors.  All  targets  are 
clearly  de'ected  by  all  the  algorithms  so  far  discussed,  and  their  spatial 

frequencies  determined  to  within  a  fraction  of  a  beamwidth.  This  behaviour  is 
typical. 

Returning  to  the  BSA  and  PIM  beamformers.  we  see  that  the  lesults  (plot 
j)  are  limited  in  application  because  of  the  high  sidelobe  levels.  Although  the 
BSA  plot.,  for  example,  contains  exactly  the  same  information  as  the  MLM 
estimate.,  it  is  not  presented  in  a  useful  fashion  for  the  location  of  multiple 
point  target  emitters.,  because  the  prior  information  that  such  sources  result  in 
a  sin(x)  /x  response  from  the  correlation  action  of  the  filter  has  not  been 
used  :  the  process  of  parameter  estimation  is  only  half  completed.  We  will 
see  that  MUSIC  and  its  related  processors  are  capable  of  effecting  this  task  of 
fitting  multiple  sin(x)/x  functions  optimally  to  the  data  when  we  examine  the 
results  obtained  from  a  focal  plane  array  which  has  an  array  manifold 
consisting  of  just  such  waveforms. 

Finally  in  this  section  we  present  the  result  of  processing  the  compacted 
data  which  results  from  the  coefficients  method  discussed  in  section  3.7.2  by 
the  MUSIC  algorithm.  The  original  data  is  exactly  the  same  as  used  for  all  the 
other  plots,  but  has  been  projected  onto  only  the  most  significant  ten 
eigenvectors  of  the  array  manifold  (corresponding  to  a  S/N  ratio  of  50  dB  in 
the  case  of  a  Shannon  number  of  six) .  The  eigen  analysis  required  for  MUSIC 
Is  thus  significantly  faster  for  the  resulting  (10  x  10)  covariance  estimate  than 
for  the  original  (20  x  20)  matrix.  Apart  fro  n  a  slight  discrepancy  in  the 
estimated  signal  powers,  the  result  plotted  in  Fig  4.5k  is  virtually 
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Fig  4  fit  Simulation  of  a  coharant  multipath  problam  Tha  thirty  snapshots  of 
data  ara  ganaratad  by  simulating  two  coharant  amitters  saparatad  Dy  0.4 
cyclas  maasurad  across  tha  twenty  sansor  apartura  Tha  rasuits  plottad  ara 
thosa  glvan  by  BSA  and  PiM 


Fig  4  6b  Tha  KTSC1  algorithm  acting  dlractly  on  tha  data  of  tha  pravious 
flgura  is  unabla  to  discrlminata  tha  two  signals. 
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Fig  4  Oc.  Utilising  a  sub-apartura  of  tan  alamants  to  pra-procasa  tha  data., 
tha  KTSC1  algorithm  succaads  In  discriminating  tha  two  coharant  signals 
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Fig  4  6d  Taking  only  a  slngla  snapshot  ot  data  from  tha  twanty  alamant 
array,  and  pra-procasslng  with  a  sub-apartura.  tha  KTSC1  algorithm  aga‘n 
succaads  in  saparatmg  tha  two  conarant  sourcas  of  Fig  4  6a  This  rasult  Is 
mors  sansltlva  to  noiss  than  that  basad  on  thirty  snapshots  of  data 

56 


A 


Indistinguishable  from  that  displayed  in  plot  a. 


,  g,  5.  PEAUNQ  WITH  COHERENT  EMITTERS 

Our  next  simulation  deals  with  perhaps  the  most  difficult  discrimination 
problem  of  all.,  that  of  two  coherent  targets  wmch  appear  to  have  a  constant 
phase  difference  at  the  antenna  phase  centre.  This  could  be  the  result  of 
coherent  multipath  reflection  from  a  fixed  surface.  We  have  simulated  e  source 
of  20  dB  and  a  reflection  of  15  dB  separated  by  0.4  conventional  Rayleigh 
beamwidths..  and  with  a  relative  phase  difference  of  rr/3  radians  measured  at 
the  centre  of  a  20  element  linoar  array  Of  half  wavelength  spaced  elements 

Fig  4.6a  shows  the  result  of  processing  twenty  snapshots  of  Information 
from  this  scenario  The  effect  of  processing  many  snapshots  in  this  way  is 
simply  to  make  the  resulting  estimate  relatively  stable  against  noise.,  and  does 
not  afford  any  improvement  in  discrimination  performance.  As  can  be  seen., 
the  BSA  and  PIM  beamscans  result  in  very  similar  responses  at  this  low  S/N 
ratio,,  and  neither  gives  a  clear  indication  of  the  presence  of  the  pair  of 
signals 

Fig  4  6b  shows  the  result  obtained  from  the  KTSC1  processor  acting  on 
precisely  the  same  data  The  sample  covariance  estimate  has  a  single  high 
eigenvalue.,  and  the  reconstruction  produces  a  single  peak  between  the  two 
sources.,  and  slightly  closer  to  the  main  signal. 

if  we  now  employ  the  sub-aperturing  technique  to  pro-process  the  data 
covariance  matrix.,  as  discussed  in  section  3.  7.  3.,  the  new  covariance  estimate 
has  two  significantly  high  eigenvalues  Forming  the  appropriate  weighted  sum 
of  noise  eigenvectors  using  the  KTSC1  algorithm  produces  the  minimum  energy 
solution  to  the  problem  given  in  Fig  4.6c.  Both  signals  are  clearly  and 
accurately  located.,  and  their  associated  powers  determined  Employing 
forward-backward  shifting  of  the  sub  aperture  gives  enhanced  noise  reiection 
and  an  even  greater  degree  of  precision  in  the  frequency  estimation,  in  fact., 
because  the  data  in  each  of  *he  snapshots  is  exactly  the  same  but  for  the 
noise  and  an  amplitude  scaling,  the  sub  aperture  may  De  applied  to  any 
individual  snapshot.,  although  with  resulting  greater  sensitivity  to  noise 
perturbations.  Fig  4.6d  shows  the  SAT  operating  successfully  on  single 
snapshot  data  from  the  even  more  difficult  situation  when  the  two  signals 
appear  to  be  in  phase  at  the  array  phase  contre. 
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Fig  4.7b  Sourest  u  U*t*d  In  Tabla  4  1.  Data  has  bssn  analytad  by  tha 
forward- backward  KTSC1  mat  hod. 
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Fig  4  7c  Sourest  at  liatsd  in  labls  4  1  Data  hat  baa n  procottsd  by 
MUSIC 
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Fig  4  7d.  Source  dataUs  are  given  In  Table  4.1.  The  resul*  Is  that  given  by 
MEM. 
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Fig  4  ?e  Details  of  emitters  are  given  In  Table  4. 1  Data  has  been  processed 
by  KTSC2. 


4.  2.  6.  MULTIPLE  TARGET  DISCRIMINATION 

Fig  4. 7  presents  the  results  of  a  number  of  simulations  which  have 
attempted  to  locate  the  powers  and  spatial  frequencies  of  seven  signal  sources 
and  two  random  phase  jamming  emissions  located  within  a  space  of  only  three 
Rayleigh  beamwidths.  using  a  twenty  element  linear  phased  array.  Processing 
has  been  carried  out  sequentially.,  to  Initially  locate  the  sources  In  space, 
followed  by  sub  aperture  analysis  of  the  resulting  time  domain  information 
Details  of  the  sources  are  given  in  Table  4  1. 

Figs  4  7a  and  4.7b  show  the  results  of  processing  by  the  KTSC1  algorithm 
on  the  forward  data  and  on  the  forward  backward  data  respectively,  in  both 
instances  the  algorithm  has  made  a  commendable  estimate  of  the  required 
parameters.  Neither  processor  has  discriminated  the  two  sources  lying  closest 
to  the  jammers  in  space  as  individual  peaks,,  although  the  subsequent  analysis 
of  the  time  series  associated  with  the  single  peak  has  accurately  determined 
the  presence  of  two  emitters  and  their  powers. 

In  this  example..  MUSIC  <  Fig  4.7c).  operating  on  an  unweighted 
combination  of  the  noise  space  basis  vectors,  has  not  performed  as  success¬ 
fully  as  KTSC1 .  It  has  completely  failed  to  register  the  presence  of  the  two 
sources  at  4  and  11  dB..  and  has  only  just  recognised  the  sources  located  at 
-1  1  cycles  mlM  (Fig  4.  7d)  displays  considerably  worse  discrimination  of  the 
sources  in  this  spatial  analysis. 

Fig  4  7e  is  Interesting,,  as  It  shows  again  the  errors  which  arise  from 
inaccurate  estimates  of  the  positions  of  high  power  sources.  This  result  has 
surprisingly  been  derived  by  the  KTSC2  processor.,  and  has  notably  failed  to 
discriminate  the  main  jamming  source.  Such  a  response  nas  resulted  in 

incorrect  matching  of  the  pseudo  Inverse  to  the  problem.,  and  the  poor 

cancellation  of  the  random  jamming  signal  in  target  directions 

Finally..  In  Fig  4.  7f,.  we  show  the  unsuitable  responses  given  by  the  BSA 
and  PIM  algorithms  to  this  data  sample  No  attempt  to  process  doppler  domain 
information  has  been  made. 

4.2.7,  FOCAL  PLANE  PROCESSING 

We  now  show  the  results  of  processing  a  'simple'  high  signal  to  noise 

problem  involving  data  received  at  seven  equally  spaced  elements  in  the  focal 

plane.  The  array  samples  spatially  at  the  Nyquist  rate  The  data  has  been 
generated  from  a  model  of  two  emitters  separated  by  0  3  beamwidths  spatially,, 
and  by  only  0  1  cycles  across  the  time  domain  aperture  in  order  to  emphasise 
the  sensitivity  of  the  successful  algorithms.  The  other  intention  here,  however, 
is  to  reveal  the  unsuitability  of  the  KTSC  algorithms  for  handling  such  data,  as 
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Fig  4. 71.  Data  ganaratad  by  a  simulation  of  tha  soureaa  llstad  In  Tabia  4  \ 
has  baan  analysed  by  PtM  and  BSA.  No  subsequent  time-series  procassing 
has  baan  attamptad. 
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Fig  4.8a  A  simulation  of  d*  a  racatvad  from  two  point  aourcas  at  0  3 
baamwidths  saparatlon  by  a  a a*sn  aiamant  wyquist  spa  cad  focai  plana  antanna 
has  baan  procassa"  by  MUSIC  Soureaa  ara  saps  rated  by  only  0  1  cyclat  In 
tha  lima  domain  apartura  Twenty  snapshou  of  data  hava  baan  usad  In  tha 
coverlanca  astimata 
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Fig  4.  8d.  Ths  data  of  Fig  4.8*  procaaaad  by  MLM. 
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Fig  4  8*  Th*  dal*  of  Fig  4  8*  pra-procaasad  by  fha  coalflclantj  maihod 
ftorm  analysis  by  KTSC1  Tha  numbar  of  spurious  paaks  1$  rsducad. 
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a  result  of  the  form  of  constraint  employed. 

Fig  4  8a  Indicates  that  MUSIC  interprets  the  data  with  ease,  despite  the 
close  proximity  of  the  sources  in  the  two  dimensional  space-time  aperture 
Spurious  peaks  are  well  suppressed  It  is  also  worth  noting  that  MUSIC  is 
effectively  interpreting  the  data  in  terms  of  sin(x)/x  functions  across  the  focal 
plane  The  response  of  the  KTSC1  processing  is  depicted  In  Fig  4  8b.  and 
displays  a  large  number  of  spurious  peaks,  in  addition  to  incorrectly  estimating 
the  locations  of  the  signals.  This  is  as  a  result  of  the  unsuitable  form  of 
constraint  employed  by  this  particular  form  of  the  signal  subspace  algorithm., 
which  requires  that  all  receiving  elements  In  the  sensor  are  equally  weighted 
and  irradiated  Both  and  MLM  (Figs  4  8c  and  4  8d)  give  superior 

results,  when  compared  with  KTSC1.  although  do  not  achieve  the  same 
accuracy  as  MUSIC  in  target  location  and  power  estimation  Again.,  neither 
algorithm,,  as  implemented  here.,  is  optimal  for  the  problem  of  focal  plane 
data 

Fig  4.  8e  shows  the  result  of  analysis  using  the  coefficients  method  of  data 
compaction  followed  by  the  KTSC1  algorithm.  Although  the  spurious  peaks  are 
pronounced.,  the  estimated  signal  powers  at  these  directions  are  relatively 
lower  than  determined  by  the  previous  analysis  of  the  data  matrix  itself.,  and 
the  locational  estimates  of  the  signal  sources  are  considerably  more  accurate 
The  coefficients  representation  of  the  data  thus  appears  to  have  made  the 
resulting  total  algorithm  less  sensitive  to  the  original  data  format  Although  it  is 
not  presented  here,  we  note  that  the  coefficients/MUSlC  method  produces 
identical  results  to  MUSIC  alone 

Finally  in  this  section.,  we  demonstrate  how  the  same  seven  element  focal 
plane  array  can  locate  at  least  seven  emitters  within  the  five  beamwidth  field  of 
view  The  scenario  consists  of  five  partially  decorrelated  sources  and  two  high 
power  random  phase  jamming  signals,,  as  listed  in  Table  4.2  Fig  4  9a  shows 
the  eigenvalue  spectrum  of  the  sample  covariance  estimate  MDL  has 
determined  a  signal  subspace  dimension  of  four.,  and  so  MUSIC  has  operated 
on  the  remaining  three  noise  subspace  eigenvectors  to  produce  the  result 
shown  in  Fig  4.9b  Four  of  the  five  directions  have  been  determined  to  within 
a  fraction  of  a  beamwidth,  and  subsequent  doppler  processing  of  the  data  has 
evaluated  the  powers  associated  with  each  of  the  signals.,  and  demonstrated 
the  excellent  rejection  of  the  lamming  signals  in  the  source  directions 
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Source 

Power 

(dB) 

Position 

(cycles) 

1 

a 

-1.5 

2 

s 

-1.5 

3 

15 

0.25 

4 

12 

0.25 

5 

25 

0.2 

6 

55 

1.2 

7 

50 

0.8 

T*di«  4.2.  Source  powara  and  locatlc.»4  uaad  lh  tha  toc»l  plan*  array 
•imuiatlon  dlacuaaad  In  auction  4.2.7  and  llluatratad  by  Fig  4 


Fig  4.9a  Tha  aoureaa  llatad  In  Tabla  4  2  hava  Daan  almu'etad  to  glva  data 
from  an  array  of  tha  iimi  apaciflcaiton  at  givan  In  Fig  4  8a  Tha  diagram 
shows  tha  alganvalua  apactrum  MDl  analyala  indicataa  that  tha  flrat  four 
aiganvactora  datma  tha  aignat  aubapaca 
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Fig  4.9b  Reconstruction  ol  the  sources  listed  In  Teble  4.2  es  attempted  by 
MUSIC  end  subsequent  time  series  enetysls. 


Fig  4  10s  The  disgrsm  shows  the  result  of  processing  a  simulation  ot  two 
5dB  sine  wave  sources,  separated  spatially  by  0  \S  cycles  and  by  0  9  cycles 
In  the  time  domain  aperture  the  twenty  snapshots  of  data  have  been  analysed 
by  the  KTSC1  method,  solving  the  equation  e,H  w  »  0.  where  e^  Is  the  matrix 
ot  spatial  domain  signal  eigenvectors  Although  the  signal  subspace  dimension 
was  clearly  two.  the  result  contains  only  a  single  peak 
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FiQ  4  10D.  Processing  the  time  series  associated  with  the  single  direction 
obtained  from  Fig  4.  10a.  using  a  moving  window/KTSCI  algorithm,  allows  ua 
to  determine  the  two  signal  frequencies. 


•  3  1.3  2.3 


Fig  4  10c.  The  data  matrix  analysed  In  Fig  4.  10a  may  alternatively  be  used  to 
determine  the  time  domain  signal  subspace  eigenvectors.  This  matrix  may 
then  be  used  to  derive  the  time  domain  estimate  shown  here  by  solving  the 
KTSC1  equation.,  xar  •  •  0.  '  * 


4.2.8.  IMPPOVING  ANGULAR  ESTIMATES  THROUGH  DOPPLER  PROCESSING 
The  results  given  In  Fig  4. 10  demonstrate  that  the  apparent  limits  of 

angular  discrimination  for  a  particular  algorithm  may  be  exceeded  the 
sources  may  first  be  discriminated  in  another  domain.,  such  as  the  doppler 
domain  In  this  example  We  have  modelled  two  5dB  sources  separated  across 
the  spatial  aperture  by  only  0  15  cycles.,  but  by  0.9  cycles  in  the  time  domain 
aperture. 

Eigenanalysis  of  the  data  covariance  estimate  gives  a  signal  subspace 
dimension  of  two  Analysis  using  the  two  singular  vectors  corresponding  to  the 
dimension  of  the  antenna  aperture.,  using  KTSC1.  produces  only  a  single 
peak,  approximately  centrally  located  between  the  two  sources  (Fig  4.  ioa) 
because  of  the  correlation  between  the  vectors  of  the  array  manifold  A  moving 
window  analysis  of  the  single  time  series  derived  from  this  directional  estimate 
reveals  the  frequencies  of  the  two  sources  with  only  a  slight  error  (Fig 
4.10b).  If.,  however.,  our  interest  is  to  locate  the  sources  spatially.,  this  result 
is  not  suitable.  We  now  show  that,,  because  of  the  decorrelation  present  in  the 
time  domain,  it  is  possible  to  find  the  directions  of  the  emissions  with  some 
degree  of  accuracy 

Instead  of  using  the  spatial  domain  singular  vectors.,  we  first  of  all  use  the 
two  signal  space  vectors  in  the  time  domain  Processing,  as  before.,  with  the 
KTSC1  algorithm,  we  determine  the  two  frequencies  to  within  a  fraction  of  a 
cycle  (Fig  4.10c).  Subsequent  analysis  of  each  of  ihe  data  series  associated 
with  these  frequencies  in  turn  enables  the  locations  of  both  of  the  emitters  to 
be  found  Such  analysis  may  be  carried  out  using  the  familiar  moving  window., 
followed  by  eigen  analysis  and  determination  of  the  appropriate  signal  subspace 
Singular  vectors  (Fig  4  lOd) .  Alternatively  in  this  situation.,  where  each 
independent  data  series  corresponds  to  only  a  single  target.,  tne  targets  may 
be  located  by  a  simple  BSA  examination  of  each  appropriate  series.  Fig  4  lOe 
shows  plots  of  the  function  l/(l-f(x))...  where  f(x)  is  the  normalised  BSA 
result.,  to  show  the  peak  positions  more  clearly.  This  latter  method  is  clearly 
much  faster  to  compute  than  the  moving  window,  but  is  obviously  not  capable 
of  finding  multiple  sources  at  a  particular  frequency  with  the  same  accuracy 

4.2.9.  ITEFtATIVE  PIM  PROCESSING 

Using  an  example  of  a  twenty  element  circular  array  of  detectors.,  looking 
only  in  the  plane  of  the  antenna.,  examples  of  the  results  of  two  different  forms 
of  iteiation  on  the  PIM  and  BSA  beamformers  are  presented  in  Figs  4.  12  and 
4. 13, 

The  target  scenario  consists  of  three  low  amplitude  sources  sufficiently 
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Fig  4. 12c.  This  figure  show*  the  first  siege  of  an  Iterated  reconstruction  of  the 
targets.  A  weighting  function  has  been  applied  to  the  original  ar'ay  manifold, 
which  assigns  unity  weight  to  those  regions  of  the  original  picture  which  fall 
above  a  threshold,  as  described  In  the  text. 
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Fig  4  12b  The  eighth  stage  of  the  Iteration  star  .ad  In  Fig  4  12t»,  Showing  the 
presence  of  a  spurious  peak  caused  by  'Incorrect*  thresholding  at  one  point 
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separated  in  aperture  that  they  are  initially  just  discriminated  by  PIM.  but  not 
quite  by  the  8SA  The  result  of  processing  the  data  by  MUSIC  using  the 
normal  pre-iteration  array  manifold  is  given  for  reference  in  Fig  4  11a..  and 
compares  with  the  equivalent  PIM  and  BSA  responses  in  Fig  4  11b.  The  latter 
figure  also  shows  the  function  1  /  <  T  — f < x>  >  derived  from  the  PIM  estimator 

Whereas  MUSIC  has  easily  determined  both  signal  locations  and  powers., 
the  BSA  algorithm  has  failed  to  produce  a  clear  Indication  of  anything.,  having 
only  a  single  broad  response  PIM  has  just  managed  to  separate  the  two 
30  dB  signals.,  and  the  possible  presence  of  the  third  slgi  jl  Is  indicated  by  a 
peak  only  slightly  above  the  level  of  the  sldelobes.,  but  at  an  incorrect 
direction. 

Fig  4  12a  shows  the  first  Iteration  of  this  scene  using  a  weighting  function 
across  the  A  matrix  consisting  of  unity  weight  for  regions  of  the  plot  which 
initially  fell  above  a  threshold  value  of  -5  dB..  and  the  same  threshold  value 
across  the  rest  of  the  picture  Estimates  of  the  locations  of  the  two  central 
emitters  are  slightly  more  accurate  than  before.,  and  the  main  sidelobe  levels 
are  slighlty  suppressed  However.,  the  unity  weighting  in  the  region  of  the  20 
dB  signal  has  not  been  applied  to  the  correct  area  because  the  corresponding 
lobe  of  the  initial  response  did  not  accurately  reflect  the  position  of  this 
source  Fig  4  12b  shows  the  result  of  the  eighth  iteration,  based  on 
intermediate  results  utilising  successively  lower  threshold  levels  We  see  that., 
although  the  plot  perhaps  looks  'cleaner'  than  before  because  of  the 
suppression  of  the  sidelobes.  no  higher  degree  of  discrimination  has  been 
achieved  In  fact.,  the  response  of  the  PIM  beamformer  has  apparently  been 
flattened  in  the  region  of  the  two  central  sources.,  and  a  spurious  peak  has 
appeared  caused  by  an  'Incorrect'  thresholding  decision  at  one  stage  of  the 
iteration  Although  the  thresholding  error  which  occured  here  was  the  result  of 
human  'error'.,.  It  could  equally  have  happened  If  the  thresholding  had  been 
carried  out  on  the  basis  of  some  other  prior  assumptions  The  only  iteration 
guaranteed  not  to  admit  such  a  possibility  of  error  is  one  supported  by  exact 
knowledge  of  the  scene  to  be  reconstructed 

Fig  4  13  shows  the  eighth  iteration  of  the  same  scene  achieved  using  a 
weighting  function  which  was  derived  from  the  square  root  of  the  amplitude  of 
each  succeeding  result  where  it  fell  above  the  threshold  level.,  and  set  to  the 
threshold  elsewhere  This  has  apparently  resulted  in  an  enhancement  of 
discrimination  performance  of  both  beamformers..  in  that  both  now  have  formed 
a  clear  estimate  of  the  directions  of  the  two  central  sources  and  the 
sidelobes  have  been  suppressed  to  very  low  levels  However.,  the  improved 
discrimination  in  the  central  region  has  primarily  occurred  because  of  the  dip 


In  the  initial  PIM  estimate  had  that  not  been  present.,  the  result  after  any 
number  of  iterations  would  have  been  only  a  single  central  peak.  We  also  note 
that  the  peak  corresponding  roughly  to  the  lower  power  target,  although  now 
clear  of  the  background.,  has  been  considerably  suppressed  itself.,  and  is  no 
more  accurately  placed  than  the  lobe  in  the  original  response  Such  a  clean 
response  has  only  been  obtained  because  the  operator  had  clear  knowledge  of 
the  locations  of  the  simulated  targets,  and  could  easily  have  resolved  spurious 
peaks  as  in  Fig  4  12b  without  such  knowledge. 

It  is  clear  that  such  iterative  processing  is  extremely  prone  to  error.,  and 
is  unable  to  cope  with  the  location  of  sources  which  faM  below  the  sidelobe  or 
interference  level  present  In  the  initial  estimate.  MUSIC.,  on  the  other  hand 
has  dealt  with  the  same  data  extremely  efficiently  and  accurately.,  on  the  basis 
of  the  assumption  that  the  targets  may  be  represented  as  point  emitters. 

For  comparison.  Fig  4.  14  shows  the  result  of  raising  the  values  obtained 
by  the  initial  PiM  estimate  to  the  power  of  eight  This  simple  operation  has 
cleariy  had  almost  as  riuch  impact  on  the  result  as  did  the  iterations  We 
conclude  that  iterative  processing  of  the  pseudo  inverse  in  this  manner  does 
not  appear  to  provide  a  useful  means  of  improving  discrimination  of 
decorrelated  point  targets 
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Fig  4.13  Tha  aighth  stags  of  an  Itsration  of  Fig  4.11b.  In  which  tho 
waightlng  function  amployad  appilss  a  walght  aqua!  to  tha  aquara  root  of  tha 
amphtuda  of  tha  pravlous  itaration  whan  this  falls  abova  tha  thrashold  lavat.  as 
dascribad  In  t  is  tart  Tha  rasult  la  a  tmoothar.  'ciaanaf  looking  raault  Tha 
BSA  rasult  has  apparantly  Improvad  discrimination,  although  this  Is  bacausa 
tha  original  PIM  was  usad  as  tha  starling  point  for  tha  procsss 
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Fig  4  14  This  diagram  shows  tha  original  BSA  rasult  of  Fig  4  11b  along  with 
tha  corraspondlng  Pim  raconstruction  which  has  basn  plotted  raisad  to  tha 
powsr  sight 


75 


S.  CONCLUSIONS 


We  have  demonstrated  the  successful  application  of  a  number  of  modern 
signal  processing  algorithms  acting  on  short  samples  of  data  from  simulated 
far  field  polnt-iike  sources  subject  to  observation  noise.  These  achieve 
improved  signal  discrimination  and  noise  reduction  over  conventional  techniques 
through  the  application  of  suitable  constraints  to  the  solution.  Such  constraints 
are  satisfied  in  a  least  squares  sense.,  which  we  have  shown  may  be  achieved 
or  understood  through  the  application  of  singular  value  decomposition. 

It  is  our  belief  that  algorithms  such  as  MUSIC  offer  worthwhile  gains  in 
both  the  quantity  and  quality  of  info-matlon  which  may  be  retrieved  from  radar 
phased  array  and  time  series  data.,  through  the  clear  and  controlled 
Introduction  of  suitable  prior  knowledge  about  the  nature  of  the  problem  to  be 
solved  We  have  demonstrated  the  algorithms  acting  on  data  of  a  fairly 
Idealistic  origin  Our  signals  have  originated  from  point  sources,  and  have 
been  subject  to  additive  Gaussian  noise.,  and  we  have  assumed  in  the  results 
presented  in  this  report  that  we  have  knowledge  of  any  perturbations  present  in 
the  array  manifolds  used  However,  the  results  obtained  strongly  suggest  that 
further  effort  should  be  directed  Into  Investigation  of  their  behaviour  when 
handling  real  data  Phased  array  radars  have  many  inherent  advantages  in  the 
form  of  robustness.,  beam  agility,,  and  so  on.,  but  it  seems  short-sighted  not  to 
seriously  experiment  with  improved  methods  of  handling  the  vast  amounts  of 
information  wh>ch  they  are  obviously  capable  of  providing  Although  the 
processing  load  is  heavy.,  we  have  demonstrated  possible  approaches  which 
might  possibly  result  in  faster  algorithms,  and  the  development  of  solid  state 
devices  is  advancing  at  such  a  rate  that  the  realisation  of  on  one  signal 
processors  performing  eigenanalysls  Is  a  real  possiblity. 

We  are  currently  evaluating  the  performance  of  these  algorithms  for  a 
variety  of  applications  such  as  resolving  multipath  duplicated  images  and 
counting  jammers  for  raid  strength  analysis.  Work  is  also  continuing  on 
possible  methods  of  improving  the  discrimination  of  sources  through  the 
inclusion  of  additional  prior  knowledge  It  is  of  vital  importance  that  the 
performance  of  such  algorithms  should  be  tested  on  digitised  signals  from  real 
phased  arrays.,  in  order  to  assess  their  sensitivity  tc  factors  which  are  not 
present  to  date  in  our  simple  simulations 

Of  the  algorithms  tested  here.  MUSIC..  EMU  and  KTSC1  appear  to  offer 
most  promise  As  we  have  shown,  they  appear  to  be  very  closely  related 
Although  by  suitable  weighting  of  the  noise  subspace  eigenvectors..  KTSC1  is 
capable  of  better  noise  rejection  than  MUSIC  when  dealing  with  short  samples 


of  data.,  we  have  demonstrated  with  the  example  of  the  focal  plane  data  in 
Section  4  2  7  that..  In  Its  current  form,,  it  Is  only  suitable  for  aperture  plane 
or  time  series  data  analysis.  MUSIC  and  the  eigenvalue  weighted  version  of  the 
algorithm,  EMU,  appear  to  be  most  widely  applicable  methods.,  sine'  they  do 
not  make  use  of  assumptions  about  the  conformation  of  the  sensor  array.  True 
maximum  likelihood  parameter  estimation  remains  a  more  distant  goal  because 
of  the  associated  weight  of  processing.  However,  it  seems  from  recently 
published  results  (Bdhme  (1985),  Hudson  (1985))  that  methods  of 
accelerated  approximate  maximum  likelihood  estimation  may  be  possible,  and  it 
is  certainly  the  case  that  such  approaches  warrant  more  detailed  Investigation 
and  characterisation. 
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